3

piモンテカルロスタイルを計算するこのランダム関数があります:

ここに画像の説明を入力してください

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 程度に制限されます。円周率の桁数を計算するには、ある種の任意精度の数値形式が必要です。

  • アルゴリズムを高速化するには、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 に答える