蒙特卡罗π计算能用于创造世界纪录吗？

36 浏览
0 Comments

蒙特卡罗π计算能用于创造世界纪录吗？

我有这个随机函数来计算圆周率蒙特卡罗风格：

\"Enter

max=10000000;
format long;
in = 0;
tic
for k=1:max
    x = rand();
    y = rand();
    if sqrt(x^2 + y^2) < 1
        in = in + 1;
    end
end
toc
calc_pi = 4*in/max 
epsilon = abs(pi - calc_pi)
calcPi(100000000);

如果我可以迭代10e100次，这个算法是否可以与世界纪录竞争？如果可以，如何找到迭代次数可以给出第N位数？

admin 更改状态以发布 2023年5月23日
0

2 答案

0 Comments

简短的回答：不行。

根据你的链接，世界记录是1e13位数字。如果您可以运行蒙特卡罗算法获得10e100个样本，那么您将获得相对RMS误差为1 / sqrt（10e100）=。3e-50的pi估计值（请参见以下内容）。这个精度仅有第50个数字的级别。此外，这只是“概率”精度：您不能确定前50个数字是正确的；您只能告诉它们很有可能是正确的。

找到蒙特卡罗样本需要多少个数字的一般规则是：M蒙特卡罗样本将为您提供1 / sqrt（M）的相对RMS精度。这意味着估计与真实值的偏差为真实值的1 / sqrt（M）部分，近似地。为了相当自信地确定第N个数字是正确的，您需要比10 ^ -N略好一些的相对RMS精度，这根据声明的规则需要M = 10 ^（2N）个样本。

因此，如果您想要（概率上的）1e13个数字的精度，则需要10 ^ 2e13个蒙特卡罗样本，这是不可管理的。

0
0 Comments

这是一个计算圆周率的不错的练习，但它可能是一个非常低效的练习。一些评论：

  • 我的统计学知识在喝咖啡之前已经生疏，但我猜错误率与 1 / sqrt(n_guess) 成比例。要获得 N 位数字，您需要一个 10^(-N) 的错误率，因此您需要大约 (10^N)^2 次随机猜测。如果你像你所提议的那样进行 1e100 次猜测，你只能获得大约 50 位数的圆周率！计算迭代次数是数字位数的某个指数函数，这非常慢。好的算法可能是您想要的数字位数的线性。

  • 由于需要大量猜测，您必须开始质疑您的随机数生成器的质量。

  • 你的算法将被浮点误差限制在 1e-16 或者这样的一个范围内。计算圆周率需要某种任意精度的数字格式。

  • 为了加快算法的速度，您可以省略 sqrt()。

  • 不要使用一个名为 max 的变量，这会覆盖一个现有的函数。使用 n_guess 等变量。

快速而肮脏的测试来证明我的理论（咖啡后）：

pie = @(n) 4 * nnz(rand(n,1).^2 + rand(n, 1).^2 < 1) / n;
ntrial = round(logspace(1, 8, 100));
pies = arrayfun(pie, ntrial);
loglog(ntrial, abs(pies - pi), '.k', ntrial, ntrial.^-.5, '--r')
xlabel('ntrials')
ylabel('epsilon')
legend('Monte Carlo', '1 / sqrt(ntrial)')

Monte Carlo result

0