3

我有这个随机函数来计算 pi Monte Carlo 风格：

在此处输入图像描述

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 位的迭代次数？

4

2 回答 2

4

这是计算 pi 的一个很好的练习，但它可能是一个非常低效的练习。一些备注：

  • 在我喝咖啡之前，我的统计数据已经生锈了，但我猜误差与1 / sqrt(n_guess). 要获得 N 位数，您需要 的误差10^(-N)，因此您需要大致(10^N)^2随机猜测。如果您按照您的建议进行 1e100 次猜测，您将只能得到 50 位 pi 的数量级！因此，迭代次数是所需位数的指数函数，这非常慢。一个好的算法可能在你想要的位数上是线性的。

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

  • 您的算法将被浮点错误限制为 1e-16 左右。计算 pi 的位数需要某种任意精度的数字格式。

  • 为了加快你的算法，你可以省略 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)')

蒙特卡洛结果

于 2013-08-09T06:43:53.533 回答
1

简短的回答：没有。

根据您的链接，世界纪录是 1e13 位数。如果您可以运行 Monte Carlo 算法来获得 10e100 个样本，您将获得 pi 的估计值，其相对 RMS 误差为 1/sqrt(10e100) = .3e-50（见下文）。这仅是第 50 位的精度。此外，它只是一种“概率”精度：您无法确定前 50 位数字是否正确；你只能说它们很有可能是正确的。

找到 N 位精度所需的蒙特卡洛样本的一般规则是：M 个蒙特卡洛样本将为您提供 1/sqrt(M) 的相对 RMS 精度。这意味着估计值与真实值的偏差大约为真实值的 1/sqrt(M) 分数。为了合理地确信第 N 个数字是正确的，您需要一个比 10^-N 稍好一点的相对 RMS 精度，根据规定的规则需要 M = 10^(2N) 个样本。

因此，如果您想要 1e13 位的（概率）精度，则需要 10^2e13 个 Monte Carlo 样本，这是无法管理的。

于 2013-08-09T07:25:16.793 回答