モンテカルロ法で円周率を求めてみるが、乱数発生に問題あり？？
モンテカルロシミュレーションを用いて円周率を求める話である。なんといっても、モンテカルロシミュレーションの命は乱数である。乱数がどれだけ不規則かで計算結果の精度が決まる。
３月１日のブログ「不規則性とは何か？ 円周率数値並びは乱数として用いるに十分に不規則か？ 」で乱数の不規則性について検討を加えた。今回は、この乱数を用いて円周率を求めてみることにした。その結果を示したのが下の図であるが、円周率よりランダムに取り出した９桁の数字、Excelで発生させた９桁の数字、そして独自の方法で作り出した９桁の数字を不規則な数値として用い、円周率を求めた。繰り返し回数を多くすると円周率の数値は円周率の真値からは約０．０００５程度離れたところまでは迫れるとの結果となった。
モンテカルロシミュレーションの場合、その計算回数を増やしていけばいくほど、得られる結果は真値に近づくはずであるが、ここが限界であった。
さて、エクセルを用いて円周率を求める話はウエブ上に多くある。表計算ソフト「Microsoft Excel」の使い方にエクセルの表画面で、モンテカルロ法を用いて円周率を求める方法が示されていたが、この方法を用いると、１００００回の繰り返しで、３．１３１２、Excel VBAも用いての１０００万回繰り返しで３．１４１５８６８となったとある。私がこのプログラムで計算した第１回目の値は3.143742、２回目の値は3.140984と、誤差０．０００５近辺まで迫っている。
なお、この文献に示されている文献には赤字で示したRandomizeがない。これがないと、常に同じ乱数列が使用されて、常に同じ結果を示すことになる。Randomizeは計算ごとに違う乱数列とするための命令語である。Randomizeがないオリジナルのプログラムで、１０００万回繰り返し時に得られた値は、文献値とは少し異なる3.1417572
であり、何度計算してもこの値となった。上の私の計算結果はプログラムにRandomizeを加えて得られた結果である。
Sub hit_miss()
' １０のｘ乗までの計算
Dim hit As Long
x = 7
Cells(1, 1) = "総試行数"
Cells(1, 2) = "円周率"
For j = 0 To x
hit = 0
For i = 1 To 10 ^ j
Randomize
If (Sqr((Rnd * 2 - 1) ^ 2 + (Rnd * 2 - 1) ^ 2)) < 1 Then
hit = hit + 1
End If
Next i
Cells(j + 2, 1) = 10 ^ j
Cells(j + 2, 2) = hit / 10 ^ j * 4
Next j
End Sub
モンテカルロシミュレーション（Wikipedia)には「30,000点をランダムに置いたとき、πの推定量は実際の値から0.07%以下の誤差の範囲内にあった。」とあり、３．１４×０．０７％＝０．００２である。乱数セットがしっかりすれば、計算精度は限りなく真値に近づくものと考えられる。現実にはエクセル発生の乱数には限界があるのだろう。ということで、本検討はここで終了とした。
モンテカルロシミュレーション（Wikipedia）より
円周率より発生させた乱数
ExcelのRnd関数で発生させた乱数
独自の方法で発生させた乱数
ブログ一覧に戻る ホームページ「アルケミストの小部屋」に戻る
３月１日のブログ「不規則性とは何か？ 円周率数値並びは乱数として用いるに十分に不規則か？ 」で乱数の不規則性について検討を加えた。今回は、この乱数を用いて円周率を求めてみることにした。その結果を示したのが下の図であるが、円周率よりランダムに取り出した９桁の数字、Excelで発生させた９桁の数字、そして独自の方法で作り出した９桁の数字を不規則な数値として用い、円周率を求めた。繰り返し回数を多くすると円周率の数値は円周率の真値からは約０．０００５程度離れたところまでは迫れるとの結果となった。
モンテカルロシミュレーションの場合、その計算回数を増やしていけばいくほど、得られる結果は真値に近づくはずであるが、ここが限界であった。
さて、エクセルを用いて円周率を求める話はウエブ上に多くある。表計算ソフト「Microsoft Excel」の使い方にエクセルの表画面で、モンテカルロ法を用いて円周率を求める方法が示されていたが、この方法を用いると、１００００回の繰り返しで、３．１３１２、Excel VBAも用いての１０００万回繰り返しで３．１４１５８６８となったとある。私がこのプログラムで計算した第１回目の値は3.143742、２回目の値は3.140984と、誤差０．０００５近辺まで迫っている。
なお、この文献に示されている文献には赤字で示したRandomizeがない。これがないと、常に同じ乱数列が使用されて、常に同じ結果を示すことになる。Randomizeは計算ごとに違う乱数列とするための命令語である。Randomizeがないオリジナルのプログラムで、１０００万回繰り返し時に得られた値は、文献値とは少し異なる3.1417572
であり、何度計算してもこの値となった。上の私の計算結果はプログラムにRandomizeを加えて得られた結果である。
Sub hit_miss()
' １０のｘ乗までの計算
Dim hit As Long
x = 7
Cells(1, 1) = "総試行数"
Cells(1, 2) = "円周率"
For j = 0 To x
hit = 0
For i = 1 To 10 ^ j
Randomize
If (Sqr((Rnd * 2 - 1) ^ 2 + (Rnd * 2 - 1) ^ 2)) < 1 Then
hit = hit + 1
End If
Next i
Cells(j + 2, 1) = 10 ^ j
Cells(j + 2, 2) = hit / 10 ^ j * 4
Next j
End Sub
モンテカルロシミュレーション（Wikipedia)には「30,000点をランダムに置いたとき、πの推定量は実際の値から0.07%以下の誤差の範囲内にあった。」とあり、３．１４×０．０７％＝０．００２である。乱数セットがしっかりすれば、計算精度は限りなく真値に近づくものと考えられる。現実にはエクセル発生の乱数には限界があるのだろう。ということで、本検討はここで終了とした。
モンテカルロシミュレーション（Wikipedia）より
円周率より発生させた乱数
ExcelのRnd関数で発生させた乱数
独自の方法で発生させた乱数
ブログ一覧に戻る ホームページ「アルケミストの小部屋」に戻る
この記事へのコメント
少し長ったらしくて恐縮ですが、
Σ（計算値－真値）／計算回数
＝Σ（計算値）／計算回数－Σ真値／計算回数
＝Σ（計算値）／計算回数－真値×計算回数／計算回数
＝Σ（計算値）/計算回数－真値
もし計算値がある値を中心としてバラツキを持っているならば、計算回数を増やすにつれて値が収束する。乱数列がしっかりとしていれば、理想的には第一項は真値に限りなく近づくことになる。
上で示した文献「表計算ソフトの使い方」では繰り返し回数と得られた値が次のようであったと記されている。
繰り返し回数 得られた円周率
１ ０
１０ ４．０
１００ ３．０８
１０００ ３．２５２
１万 ３．１２２８
１０万 ３．１４４４０
１００万 ３．１４２０９６
１０００万 ３．１４１５８６８