1 # モ ン テ カ ル ロ 法 で 円 周 率 を 求 め る 2 i m p o r t r a n d o m
3 r a n d o m . s e e d (0) 4
5 def c o u n t ( n ):
6 " n個 の 点 を[ -1 ,1]*[ -1 ,1]に ば ら ま い て ， 円 の 中 に あ る 点 の 個 数 を 数 え る 関 数"
7 p =0 # 以 下 で カ ウ ン ト す る 値 を pに 入 れ る 8 for i in r a n g e( n ):
9 x = r a n d o m . u n i f o r m ( -1 ,1) # 乱 数 の x座 標 10 y = r a n d o m . u n i f o r m ( -1 ,1) # 乱 数 の y座 標
11 if x ^2+ y ^2 <1: # も しx ^2+ y ^2 <1な ら ば
12 p = p +1 # pを 一 つ 増 や す
13 r e t u r n p # カ ウ ン ト し た 点 の 数 pを 返 す
14
15 n = 1 0 0 # サ ン プ ル 点 の 個 数
16
17 pai = 4 . 0 * c o u n t ( n )/ n # こ れ が 円 周 率
18 p r i n t pai # p a iを 出 力 す る
実行結果の例 3.28000000000000
上のプログラムを実行してみましょう。出てきた結果が円周率の近似値です。サンプル点nの値を1から
10000ぐらいの範囲で変化させて，どのぐらい円周率の真の値3.1415926535· · · に近くなるか試してみましょ
う。やってみると上のプログラムはそれほど正確な円周率の値を出さないことが分かります。
次に，上のプログラムでnについて計算した円周率をpai(n)と定義して，これをnを変化させた結果を見 てみましょう。
1 # こ こ ま で は 上 の プ ロ グ ラ ム 2 # の1 3行 目 ま で と 同 じ に す る 。 3
4 def pai ( n ):
5 r e t u r n 4 . 0 * c o u n t ( n )/ n 6
7 for i in r a n g e(1 ,51):
8 p r i n t 4 0 0 * i , pai ( 4 0 0 * i )
実行結果の例 400 3.07000000000000
800 3.05000000000000 1200 3.14000000000000
……
19200 3.13895833333333 19600 3.11795918367347 20000 3.13840000000000
20000個のサンプル点を取って計算してもそれほ0.01の桁が異なっています。これは確率論の基本的な事実
で，分散を計算すればすぐにわかります。この方法ではサンプル点の数をnとしたときに，円周率の誤差はお よそ1/√
nの割合になります。
さて，上のプログラムを改良して収束の様子をグラフで表してみましょう：
1 lis = []
2 for i in r a n g e(1 , 2 0 0 ) :
3 lis . a p p e n d ([ i , pai ( 1 0 0 * i )]) # ペ ア[ i , pi ( 1 0 0 * i )]の リ ス ト を 作 る 4 p1 = l i s t _ p l o t ( lis , p l o t j o i n e d =True, l e g e n d _ l a b e l =’ m o n t e c a r l o ’)
50.2 モンテカルロ法 50 乱数
5 p2 = l i n e ([(0 , pi ) ,(200 , pi )] , c o l o r =’ red ’) # 円 周 率 の 真 の 値 6 ( p1 + p2 ). s h o w ()
計算には少し時間がかかりすぎる場合は上の100*iを10*iに変えてから実行してみましょう。次のような グラフが表示されます：
0 50 100 150 200
2.9 2.95 3 3.05 3.1 3.15 3.2
3.25 monte carlo
50.2.2 収束の早さの比較
さて円周率は次のような交代級数で計算することもできます：
π= 4 arctan(1) = 4
∫ 1 0
1
1 +x2dx= 4
∫ 1 0
∑∞ n=0
(−1)nx2ndx= 4
∑∞ n=0
(−1)n 1 2n+ 1
= 4 (
1−1 3 +1
5 −1 7+· · ·
)
この級数のn番目までの和の円周率の誤差は明らかに1/n程度です。このアルゴリズムを用いて計算して円 周率と先ほどの乱数で計算した円周率との収束性を比較してみましょう：上のプログラムに続いて次を入力し ます：
1 # a r c t a nの 展 開 級 数 で 求 め た 円 周 率 の 収 束 の 様 子 を み る 2 def p i _ a r c t a n ( n ):
3 " a r c t a n (1)の 級 数 展 開 の n個 の 和 に よ り 近 似 し た 円 周 率"
4 at a n =0
5 for j in r a n g e( n ):
6 at a n = a t a n + ( 1 . 0 / ( 2 * j + 1 ) ) * ( ( - 1 ) ^ j ) 7 r e t u r n 4* a t a n
8
9 l i s 2 =[]
10 for i in r a n g e( 2 0 0 ) :
11 li s 2 . a p p e n d (( i , p i _ a r c t a n (5* i ))) 12
13 p3 = l i s t _ p l o t ( lis2 , p l o t j o i n e d =True, c o l o r =’ g r e e n ’, y m i n =+3 , y m a x = 3 . 3 ) 14 ( p1 + p2 + p3 ). s h o w ()
0 50 100 150 200 3
3.05 3.1 3.15 3.2 3.25
3.3 monte carloseries of arctan
arctan(1)の級数展開で計算した円周率の収束も必ずしも良いとは言えませんが，計算すればするほど，真
の値に近づく事が分かります。一方モンテカルロ法は収束の速度は遅く，なかなか真の値に近づかない事が分 かります。モンテカルロ法で誤差∆の割合で結果を得るにはおよそ1/∆2個のサンプル点を取る必要があり ます。モンテカルロ法が効力を発揮するのは高次元の積分値を求めるときです。
50.3 4 次元球の体積のモンテカルロ法による計算
4次元球とはR4の集合
B4={x= (x1, x2, x3, x4)∈R4|x21+x22+x23+x24≤1} (34) です。ここでは[−1,1]4の中に点をばらまいて，B4の中にある点を数える事によってB4の体積を計算し ます：
(ばらまいた点の数) : (B4の中にある点の数)∼= 24:B4の体積 なので
B4の体積∼= 16×B4の中にある点の数 ばらまいた点の数
です。したがって4次元球の体積を計算するためには次のようなプログラムを書けばよいでしょう：
1 ’ モ ン テ カ ル ロ 法 で4次 元 単 位 球 の 体 積 を 求 め る’ 2 i m p o r t r a n d o m ; r a n d o m . s e e d (0)
3 def c o u n t ( n ):
4 p =0 # カ ウ ン ト す る 数 を pと す る 5 for i in r a n g e( n ):
6 x = r a n d o m . u n i f o r m ( -1 ,1) # 乱 数 の x座 標 7 y = r a n d o m . u n i f o r m ( -1 ,1) # 乱 数 の y座 標 8 z = r a n d o m . u n i f o r m ( -1 ,1) # 乱 数 の z座 標 9 w = r a n d o m . u n i f o r m ( -1 ,1) # 乱 数 の w座 標
10 if x ^2+ y ^2+ z ^2+ w ^2 < 1: # も しx ^2+ y ^2+ z ^2+ w ^2 <1な ら ば
11 p = p +1 # pを 一 つ 増 や す