 |
本章では、モンテカルロ法とよばれる、乱数を用いて問題を解く方法を学びます。解析的に計算が困難な問題でも、制限なく解くことができます。１桁精度をあげるために、それまで要した計算の１００倍の時間がかかり、実行時間が遅いのが欠点です。
平行線と針の交わる確率
１８世紀のフランスの自然科学者 Comte de
Buffonは、πの近似値を次のような実験で求めました。まず、平面上に平行線を何本か引き、その間隔と同じ長さの針をその上に落とします。針はランダ
ムに落ち、平行線と交わるものと、そうでないものがあります。針の交わる確率が、2/πとなることから、落とした針の本数と、交わった針の本数の比とか
ら、πを求めようとするものです。
| 図５．１．１ ビュフォンの針 |
 |
針の交差する理論的な確率は、図５．１．１から次のように計算できます。落ちた針の中央の座標を x0, y0 とし、針が x 軸となす角度を θ とします。交差する確率は y 軸とは無関係で、 x0 と θ だけの関数となります。さらに、 x0 が 0 と 1/2 の間、 θ が 0 から π/2 の値の間だけを考えれば、あとはその繰り返しになります。
平行線と交わるための条件は、図５．１．１ａから x0≦cosθ となりますから、図５．１．１ｂの cosθ の下側であれば交わることとなります。この斜線部分の面積（ =1/2 ）と長方形の面積（ =π/4 ）との比が、交差する確率となります。すなわち
交差する確率 =2/π
となります。何本も針を落として見て、交わっている針の数をかぞえ、その確率がもとめられたとします。その確率からこの式に従ってπを求めようとするものです。
シミュレーションで再現するビュフォンの針
プログラムS0510.javaは、シミュレーションで針を落とす実験を再現したもので、ディスプレイ画面の上に９本の平行線を引き、そ
こに針を落とす様子が次々に表示されます。プログラムの中の（１）は針の先端の x
座標の整数部が異なっていることにより、平行線と交わっていることを検出しています。
実行してみましょう。この結果、π の値は3.14に近い数値が得られますか。（図５．１．２）
| 図５．１．２ ビュフォンの針 実行結果 |
 |
| S0510.java ビュフォンの針 | Download | 実行 | |
/* S0510.java
* ビュフォンの針
* (C) H.Ishikawa 2008
*/
package simulation;
import java.applet.*;
import java.awt.*;
import java.awt.event.*;
import window.Window;
public class S0510 extends Applet implements ActionListener {
Button button0;
Button button1;
boolean sw = false;
public void init() {
button0 = new Button(" 実行 ");
button1 = new Button("クリア");
add(button0);
add(button1);
button0.addActionListener(this);
button1.addActionListener(this);
}
public void actionPerformed(ActionEvent e) {
String label = e.getActionCommand();
if (label.equals(" 実行 "))
sw = true;
else
sw = false;
repaint();
}
public void paint(Graphics g) {
Window w;
w = new Window();
int height;
int width;
int NUMBER = 1000; /* 落とした針の本数 */
double MAX = 8.0;
double P = 0.5; /* 確率p */
int j; /* forのカウンタ */
int xl; /* forのカウンタ */
int n1 = 0; /* 針が交わった回数 */
double x0, y0; /* 針の中央の座標 */
double x1, y1, x2, y2; /* 針の先端の座標 */
double th; /* θ */
/* スクリーンの大きさを読み取る */
width = getSize().width;
height = getSize().height;
/* グラフィックの準備 */
w.setWindow(0, 0, 0, MAX, MAX * height / width,
0, height, width - 1, 0);
for (xl = 0; xl <= MAX; xl++) {
w.line(0, xl, 0, xl, MAX, g); /* 平行線 */
}
g.setColor(Color.green); /* 緑 */
/* メイン */
if (sw == true) {
for (j = 1; j <= NUMBER; j++) {
x0 = MAX * Math.random();
y0 = MAX * height / width * Math.random();
th = 2 * Math.PI * Math.random();
x1 = x0 + Math.cos(th) / 2;
y1 = y0 + Math.sin(th) / 2;
x2 = x0 - Math.cos(th) / 2;
y2 = y0 - Math.sin(th) / 2;
w.line(0, x1, y1, x2, y2, g);
if ((int)x1 != (int)x2) {
/* (1)針の両端の座標の整数部が異なっている時 */
n1 = n1 + 1;
}
}
g.setColor(Color.black);
g.drawString("pi = " + (double)(j) /(double)(n1) * 2.0,
20, 10);
}
}
}
|
大数の法則の応用
ビュフォンの針は、特別な方法でπを求めましたが、もっとこれを一般化してみます。３．４節の
大数の法則で述べたように、コイン投げを多数回（ m 回）試みると、成功の回数（ k
）との比は一定値に近づきます。モンテカルロ法は、コインを投げる代わりにある実験を行い、その成功の数を数える実験的な数学です。 π
を求める計算をシミュレーション実験し、モンテカルロ法を理解しましょう。
正方形の上に雨が降る
図５．２．１のように、半径１の１／４円を、一辺が１の正方形に接するように書くとしましょう。この正方形を、雨の降る中にさらしたとします。雨滴は均
一に、ランダムに降りそそぎますから、１／４円の内側（図のＡの部分）も、残りの部分（図のＢの部分）もぬらします。Ａの部分に落ちる雨滴の数と、正方形
の雨滴の数の比は、大数の法則により面積の比に近づくはずです。すなわち、
 |
・・・・・ |
式5.2.1 |
となります。ここで、 右辺は、π／４ですから、もし雨滴の数を数えることができれば、
 |
・・・・・ |
式5.5.2 |
として、πが計算できるはずです。
| 図5.2.1 １／４円 |
 |
プログラムのポイント
Math.random()メソッドより、区間 [0, 1) の乱数を２個発生し、それぞれを x0, y0 とします。もし
 |
・・・・・ |
式5.2.3 |
ならば、Ａの部分の雨滴に相当します。多数（ j 個）の乱数を発生させ、Ａにある回数 k を数え、
 |
・・・・・ |
式5.2.4 |
としてπを求めます。
プログラムS0520.javaは、雨滴が正方形をぬらしていくありさまをディスプレイするようにしたものです。図５．２．２のように
j が大きくなるにつれて、 π に収束していくことがわかります。
| 図５．２．２ モンテカルロ法によるπの計算結果 |
 |
| S0520.java モンテカルロ法によるπの計算 | Download | 実行 | |
/* S0520.java
* モンテカルロ法によるπ
* (C) H.Ishikawa 2008
*/
package simulation;
import java.applet.*;
import java.awt.*;
import java.awt.event.*;
import window.Window;
public class S0520 extends Applet implements ActionListener {
Button button0;
Button button1;
boolean sw = false;
public void init() {
button0 = new Button(" 実行 ");
button1 = new Button("クリア");
add(button0);
add(button1);
button0.addActionListener(this);
button1.addActionListener(this);
}
public void actionPerformed(ActionEvent e) {
String label = e.getActionCommand();
if (label.equals(" 実行 "))
sw = true;
else
sw = false;
repaint();
}
public void paint(Graphics g) {
Window w;
w = new Window();
int height;
int width;
int NUMBER = 10000; /* 雨滴の数 */
int SPACE = 30;
int j; /* forのカウンタ */
int xl; /* forのカウンタ */
int n1 = 0; /* 円の中に入った雨滴 */
double x0, y0; /* 雨滴の座標 */
double pi = 0.0; /* πの近似値 */
/* スクリーンの大きさを読み取る */
width = getSize().width;
height = getSize().height;
/* グラフィックの準備 */
w.setWindow(0, 0, 0, 1, 1, SPACE,
height/2 + SPACE, height/2 + SPACE, SPACE);
w.axis(0, "", 1.0,"", 1.0, g);
w.setWindow(1, 0, 0, NUMBER, 4,
SPACE, height - SPACE, width - SPACE, height/2 + 2*SPACE);
w.axis(1, "j", NUMBER / 10, "pi",1.0, g);
/* メイン */
if (sw == true) {
w.moveTo(1, 0.0, 0.0, g);
for (j = 1; j <= NUMBER; j ++) {
x0 = Math.random();
y0 = Math.random();
if (x0 * x0 + y0 * y0 < 1.0) {
/* 円の中に入った場合 */
n1 = n1 + 1;
g.setColor(Color.blue);
w.putPixel(0, x0, y0, g);
}
pi = (double)n1 / j * 4;
g.setColor(Color.green);
w.lineTo(1, (double)j, pi, g);
}
g.setColor(Color.black);
g.drawString(" pi = " + pi, 10, 10);
stop();
}
}
}
|
球の体積
モンテカルロ法による数値積分は、多重積分において、特に威力を発揮します。球の体積をモンテカルロ法により求めてみましょう。 半径１の球は、図５．３．１のように、
 |
・・・・・ |
式5.3.1 |
で表わされます。この球の体積 V は、図５．３．１の x-y 平面、 y-z 平面、 z-x 平面で切りとられる体積の８倍ですから、
 |
・・・・・ |
式5.3.2 |
という定積分を行うこととなります。
| 図５．３．１ 球 |
 |
プログラムのポイント
５．２の場合と同様に考え、Math.random()メソッドにより、区間 [0, 1) の一様乱数を３個発生させ、それぞれを、 x0 , y0 , z0 とします。もし、
x0 * x0 + y0 * y0 + z0 * z0 < 1
であれば、球の内側の乱数であったので、その数を数えます。全乱数の数 j と、内側であった数 k との比から、体積 V は、
v = 8 * k / j
として求めます。
プログラムS0530.javaは球の体質を、モンテカルロ法により求めるプログラムです。実行し、理論値４／３・π（約4.19）と比較しましょう。
| 図５．３．２ モンテカルロ法による球の体積計算結果 |
 |
| S0530.java モンテカルロ法による球の体積計算 | Download | 実行 | |
/* S0530.java
* モンテカルロ法による球の体積計算
* (C) H.Ishikawa 2008
*/
package simulation;
import java.applet.*;
import java.awt.*;
import java.awt.event.*;
import window.Window;
public class S0530 extends Applet implements ActionListener {
Button button0;
public void init() {
button0 = new Button(" 再実行 ");
add(button0);
button0.addActionListener(this);
}
public void actionPerformed(ActionEvent e) {
String label = e.getActionCommand();
repaint();
}
public void paint(Graphics g){
Window w ;
w = new Window();
int SPACE = 20;
int HIGHT = 400;
int WIDTH = 640;
int NUMBER = 1000; /* 乱数の組の数 */
int j; /* forのカウンタ */
int n1 = 0; /* 球の中に入った数 */
double x0, y0, z0; /* 点の座標 */
double v = 0.0; /* 体積の近似値 */
w.setWindow(0, 0, 0, NUMBER, 8,
SPACE, HIGHT-SPACE, WIDTH-SPACE, SPACE);
w.axis(0, "j", NUMBER / 10, "v", 1, g);
g.setColor(Color.green);
w.moveTo(0, 0.0, 0.0, g);
for (j = 1; j <= NUMBER; j ++) {
x0 = Math.random();
y0 = Math.random();
z0 = Math.random();
if (x0 * x0 + y0 * y0 + z0 * z0 < 1.0) {
/* 球の中に入った場合 */
n1 = n1 + 1;
}
v = (double)n1 / j * 8;
w.lineTo(0, j, v, g);
}
g.setColor(Color.black);
g.drawString(" v = " + v, 10, 10);
stop();
}
}
|
１． モンテカルロ法により、
を求めなさい（理論値 =1 - e-1 = 0.632121 ）。
（解答） A0510.java | 実行 |
２．
という定積分を行うとき、区間 [0, 1) の乱数 ri を用い、それを f(x) に次々代入し、
としても求められます。この方法により、１．と同じ定積分をもとめなさい。
（解答） A0520.java | 実行 |
３． ２．を拡張して
という多重積分を求めることができます。すなわち、３個１組の区間 [0, 1) の乱数 x0, y0, z0 を n 組発生させ、
として求めます。この方法により、
を求めなさい。
（解答） A0530.java | 実行 |
|
 |
|
