このページをはてなブックマークに追加このページを含むはてなブックマーク このページをlivedoor クリップに追加このページを含むlivedoor クリップ

はじめに

モンテカルロ法

 紙の上に適当な図形を描く。次に粒の揃った砂を用意して、紙の上に均等にばらまく。その後、図形の中と外の砂粒の数を数える。紙全体の面積に図形の中に入った砂粒の数の割合を掛ければ、その図形のおおよその面積を求めることができる。こうしたアルゴリズムをモンテカルロ法という。

 計算が難しい問題を乱数を利用して解くというアイデアは、第2次世界大戦中のアメリカのロスアラモス研究所で、原爆の研究の中で発見された。サイコロを振って答えるような方法なので、カジノで有名なモナコのモンテカルロにちなんで名付けられた。

モンテカルロ法を利用して円の面積を解く

 次のようにNとnを定義する。

  • 全体の点をN個
  • 1/4円内の点n個

 このとき、四角形の面積はr^2、1/4円の面積は\frac{n}{N}~r^2である。

 ランダム値を生成し、(x,y)の点を作り出す。それが円の中にあるかどうかを判断して、点の数を数える。円の中に入っているかどうかは三平方の定理を使う。点(x,y)がx^2+y^2~\le~r^2を満たせば、円の中ということになる。そして、乱数は0以上の数なので、x,yは0以上になる範囲で、最大値が円の半径になるような四角形の領域に点を打って考える。こうして求めた面積は1/4円なので、後で4倍すればよい。

 また、整数の乱数だと、半径の小さな円では精度が悪くなるので、double型の乱数を利用した。

VB.NETによる実験

 モンテカルロ法による円の面積を求める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
Public Class Form1
    Private Const randomNum As Integer = 100000 'ばらまく乱数の個数
 
    Private Sub Button1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button1.Click
        Dim inNum As Integer = 0  '1/4円に含まれた座標の数
        Dim x, y, pai As Single
 
        Dim width As Integer = PictureBox1.Width
        Dim height As Integer = PictureBox1.Height
        Dim g As Graphics = PictureBox1.CreateGraphics()
        g.TranslateTransform(0, height)
        g.ScaleTransform(1, -1)
 
        '--------------------------------------------------------------------------------------
        '(0,0),(1,0),(0,1),(1,1)に囲まれた領域上の座標をランダムに選択する。
        '--------------------------------------------------------------------------------------
        Randomize() '擬似乱数生成の種を生成する。
        Dim i As Integer
        For i = 0 To randomNum
            x = Rnd()
            y = Rnd()
 
            'ランダムに生成した座標が1/4円に含まれているかどうかを調べる。
            If x * x + y * y < 1 Then
                g.FillRectangle(Brushes.Blue, x * width, y * height, 1, 1)    '点を青色でべた塗りする。
                inNum += 1
            Else
                g.FillRectangle(Brushes.White, x * width, y * height, 1, 1)    '点を白色でべた塗りする。
            End If
        Next
 
        '--------------------------------------------------------------------------------------
        'πを計算する。
        '--------------------------------------------------------------------------------------
        pai = 4 * inNum / randomNum   'π/4 : 1 = inNUm : randomNumより
        ListBox1.Items.Add(pai)
    End Sub
End Class

 実行結果は次の通りである。

 10万個の点では小数点第2位までの3.14ぐらいしか一致しないことがわかる。

Javaによる実験

 モンテカルロ法による円の面積を求めるJavaのプログラム(MonteCarlo.java)を作ると、次のようになる。

Everything is expanded.Everything is shortened.
  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
 
 
 
 
 
 
 
-
|
|
|
|
-
|
|
-
|
|
|
|
|
-
|
|
!
!
|
|
|
|
!
|
|
-
|
|
|
|
|
|
-
-
|
|
-
|
|
|
!
!
|
|
|
|
|
|
|
|
|
|
|
!
!
/*
 * モンテカルロ法による円の面積計算
 * 入力:半径、試行回数
 * 出力:円の面積
 */
import java.util.*;
 
public class MonteCarlo {
    //モンテカルロ法で円の面積を求める
    //r:円の半径
    //tno:点を打つ回数
    //戻り値:面積
    static double MonteSquare(double r,int tno){
        Random rd=new Random();
        int cnt=0;
        for(int i=0;i<tno;i++){
            //x,y座標をランダムに求める
            double x=rd.nextDouble()*r;
            double y=rd.nextDouble()*r;
            
            //三平方の定理で円の中かどうかを判定
            if(x*x+y*y<=r*r){
                //円の中ならカウント
                cnt ++;
            }
        }
        
        //点の数から面積を求める(4倍する)
        double m=4*r*r*(double)cnt/(double)tno;
        return m;
    }
    
    //main
    public static void main(String[] args) {
        //乱数を発生させる回数(点の個数)
        int tryalNo=1000;
        //円の半径
        double radius=10;
        
        //引数のチェック
        if(args.length==2){
            try{
                radius=Integer.parseInt(args[0]);
                tryalNo=Integer.parseInt(args[1]);
            }catch(NumberFormatException e){
                System.out.println("モンテカルロ法で円の面積を求める");
                System.out.println("MonteCarlo <r> <n> r:円の半径 n:試行回数");
                return;
            }
        }
        
        //モンテカルロ法で円の面積を求める
        double mont=MonteSquare(radius,tryalNo);
        
        //計算で求めた値とモンテカルロ法で求めた値を比較して、表示する
        double circle=radius*radius*Math.PI;
        System.out.println("計算値(π×r×r)="+circle);
        System.out.println("モンテカルロ="+mont);
        double ratio=(1.0-mont/circle)*100;
        String rs=new java.text.DecimalFormat("0.0####").format(ratio);
        System.out.println("誤差(1-b/a):"+rs+"%");
    }
}

 クラスメソッド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)は次のようになる。

Everything is expanded.Everything is shortened.
  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
 
 
 
 
 
 
 
-
|
|
|
|
-
|
|
-
|
|
|
|
|
-
|
|
!
!
|
|
|
|
!
|
|
-
|
|
|
|
|
|
-
-
|
|
-
|
|
|
!
!
|
|
|
|
-
|
!
|
|
|
|
|
|
|
|
|
!
!
/*
 * モンテカルロ法による円の面積計算
 * 入力:半径、試行回数
 * 出力:円の面積
 */
import java.util.*;
 
public class MonteCarlo {
    //モンテカルロ法で円の面積を求める
    //r:円の半径
    //tno:点を打つ回数
    //戻り値:面積
    static double MonteSquare(double r,int tno){
        Random rd=new Random();
        int cnt=0;
        for(int i=0;i<tno;i++){
            //x,y座標をランダムに求める
            double x=rd.nextDouble()*r;
            double y=rd.nextDouble()*r;
            
            //三平方の定理で円の中かどうかを判定
            if(x*x+y*y<=r*r){
                //円の中ならカウント
                cnt ++;
            }
        }
        
        //点の数から面積を求める(4倍する)
        double m=4*r*r*(double)cnt/(double)tno;
        return m;
    }
    
    //main
    public static void main(String[] args) {
        //乱数を発生させる回数(点の個数)
        int tryalNo=1000;
        //円の半径
        double radius=10;
        
        //引数のチェック
        if(args.length==2){
            try{
                radius=Integer.parseInt(args[0]);
                tryalNo=Integer.parseInt(args[1]);
            }catch(NumberFormatException e){
                System.out.println("モンテカルロ法で円の面積を求める");
                System.out.println("MonteCarlo <r> <n> r:円の半径 n:試行回数");
                return;
            }
        }
        
        //モンテカルロ法で円の面積を求める(ここを修正)
        double mont=0;
        int rep=10;    //繰り返し数
        for(int i=0;i<rep;i++){
            mont += MonteSquare(radius,tryalNo);
        }
        mont /= rep;    //平均を求める
        
        //計算で求めた値とモンテカルロ法で求めた値を比較して、表示する
        double circle=radius*radius*Math.PI;
        System.out.println("計算値(π×r×r)="+circle);
        System.out.println("モンテカルロ="+mont);
        double ratio=(1.0-mont/circle)*100;
        String rs=new java.text.DecimalFormat("0.0####").format(ratio);
        System.out.println("誤差(1-b/a):"+rs+"%");
    }
}

参考文献

  • 『Eclipseによる体験学習 Javaではじめるアルゴリズム入門』