紙の上に適当な図形を描く。次に粒の揃った砂を用意して、紙の上に均等にばらまく。その後、図形の中と外の砂粒の数を数える。紙全体の面積に図形の中に入った砂粒の数の割合を掛ければ、その図形のおおよその面積を求めることができる。こうしたアルゴリズムをモンテカルロ法という。
計算が難しい問題を乱数を利用して解くというアイデアは、第2次世界大戦中のアメリカのロスアラモス研究所で、原爆の研究の中で発見された。サイコロを振って答えるような方法なので、カジノで有名なモナコのモンテカルロにちなんで名付けられた。
次のようにNとnを定義する。
このとき、四角形の面積は、1/4円の面積は
である。
ランダム値を生成し、(x,y)の点を作り出す。それが円の中にあるかどうかを判断して、点の数を数える。円の中に入っているかどうかは三平方の定理を使う。点(x,y)がを満たせば、円の中ということになる。そして、乱数は0以上の数なので、x,yは0以上になる範囲で、最大値が円の半径になるような四角形の領域に点を打って考える。こうして求めた面積は1/4円なので、後で4倍すればよい。
また、整数の乱数だと、半径の小さな円では精度が悪くなるので、double型の乱数を利用した。
モンテカルロ法による円の面積を求めるVB2008のプログラムを作ると、次のようになる。Form1にはListBox,PictureBox,Buttonをそれぞれ1つずつ貼り付ける。このプログラムは10万個の点をランダムに生成して、モンテカルロ法を用いてπを計算している。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 |
|
実行結果は次の通りである。
10万個の点では小数点第2位までの3.14ぐらいしか一致しないことがわかる。
モンテカルロ法による円の面積を求めるJavaのプログラム(MonteCarlo.java)を作ると、次のようになる。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
| - | | | | - | | - | | | | | - | | ! ! | | | | ! | | - | | | | | | - - | | - | | | ! ! | | | | | | | | | | | ! ! |
|
クラスメソッドMonteSquareの内部では、tnoとcntが使われている。tnoは点を打つ回数、cntは点の数を数える整数である。int型で割り算をすると、小数点以下が切り捨てられてしまうので、double型に型変換する。そうしてから円の面積を求める。
このプログラムには、コマンドラインの引数として円の半径と試行回数を指定できる。指定がなければ、半径は10、試行回数は1,000で計算される。最後に、計算で求めた円の面積と比較して、誤差をパーセントで表示するようにしてある。この誤差の値が小さい方が、より正確ということを意味する。
例:MonteCarlo.javaの実行結果
>java MonteCarlo 10 1000000 計算値(π×r×r)=314.1592653589793 モンテカルロ=314.4392 誤差(1-b/a):-0.08911%
>java MonteCarlo 10 1000000 計算値(π×r×r)=314.1592653589793 モンテカルロ=314.1732 誤差(1-b/a):-0.00444%
>java MonteCarlo 10 1000000 計算値(π×r×r)=314.1592653589793 モンテカルロ=314.012 誤差(1-b/a):0.04688%
>java MonteCarlo 10 1000000 計算値(π×r×r)=314.1592653589793 モンテカルロ=314.0484 誤差(1-b/a):0.03529%
第2引数が大きいほど、打つ点が多いので、誤差が小さくなる。また、それぞれの結果の誤差が結構ばらつきがあるので、10回分のモンテカルロ法で解いた結果の平均を考えることで、値が安定してくる。ただし、計算回数が多くなるので、その分出力を得る時間が重くなる。そのために、MonteSquareメソッドを実行している行に、次のコードを追加すればよい。
double mont=0; int rep=10; //繰り返し数 for(int i=0;i<rep;i++){ mont += MonteSquare(radius,tryalNo); } mont /= rep; //平均を求める
修正後のプログラム(MonteCarlo2.java)は次のようになる。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
| - | | | | - | | - | | | | | - | | ! ! | | | | ! | | - | | | | | | - - | | - | | | ! ! | | | | - | ! | | | | | | | | | ! ! |
|