モンテカルロ法入門

モンテカルロは カジノ で有名なモナコの首都で，その名前をとった数値計算の技法が モンテカルロ法である。つまり，いい加減に定義すると，賭博の原理を使った計算方法というわけだ．賭博というものに手を出すと一生に不幸になるので，大金持ちやその筋の人でない諸君は手を出さないように！

積分を評価する

計算の例—円周率 $\pi$ を求める

計算の原理

原点に座標を持つ半径1の円を考える。周知のように この面積は $\pi/4$ である。一方， この円の方程式は $x^2+y^2=1$ であるから， この面積は， 関数$y = \sqrt{1-x^2}$の区間 $[0,1]$の範囲と$x$軸の間に挟まれた部分の 面積になっている。 今，$\pi$ が未知として，この値をモンテカルロ法で求めよう。

上の図で $0\ge x \ge 1 $ と $0\ge y \ge 1 $ に挟まれた正方形の領域の面積は 1 である。この領域に 入るように乱数で点の位置を決めて 放り込んでいくことにする。 放り込まれた乱数 の一部は，青で示された円周の内側の領域に入ることになり，その確率は 面積の比から $1/4 \pi$ になる。

そこで，$n$個の乱数による点を正方形の領域に放り込んだとすると， 青い領域に入る点の数 $m$ は次のようになることが期待される。

$m = 1/4 n \pi$

これを変形して次の式が得られる。

$\pi = \frac{4m}{n}$

つまり，$n$ 個の 点をランダムに正方形の領域に打ち込んで， 円周の内側に入ったものをカウントして $m$ とすると， 上の式を使って $\pi$ の近似値が求められることになる。

ここで，円周の内側に入ったものだけをカウントするには， 乱数で生成された座標を $(x,y)$ として， 次の不等式を満たすものを数えていけばよい。

$ y \lt \sqrt{1-x^2} $

Ruby で書いたプログラム

次のソースは Ruby で上記のモンテカルロ法による $\pi$ の計算を 行うためのものである。走らせるには，プログラム名の後にコマンドライン引数で 試行回数を与える。また，他の言語との比較のために，計算にかかった時間が 表示される。

#!/usr/local/bin/ruby
# calcPi.rb
include Math
hit = 0
n_trial = ARGV[0].to_i
time1 = Time.now
n_trial.times{
  x = rand()
  y = rand()
  hit += 1 if y < sqrt(1-x**2)
}
time2 = Time.now
printf("%12.9f\n",4.0*hit/n_trial)
printf("%8.3f sec.",(time2-time1))

実行には次のようにコマンドラインから入力する。

>  ./calcPi.rb 1000000 (回数は適宜増減する)

Java で書いたプログラム

Java はコンパイラ言語なので，実行は Ruby よりも格段に速い。そのソースを 以下に示す。

import java.util.Random;
import java.lang.*;
public class calcPi{
  public static void main(String s[]) {
      long i,n_trial,hit=0;
      double x,y;
      n_trial = Integer.parseInt(s[0]);
      Random rand = new Random(10);

      long time1, time2;
      time1 = System.currentTimeMillis();

      for(i=0;i < n_trial;i++){
	  x = rand.nextDouble();
	  y = rand.nextDouble();
	  if(y < Math.sqrt(1-x*x))
	      hit ++;
      }
      time2 = System.currentTimeMillis();
      System.out.printf("%12.9f\n",4.0*hit/(double)n_trial);
      System.out.printf("%8.3f sec.\n",((time2-time1))/1000.0);
  }
}

コンパイルして実行するには，次のようにコマンドラインから入力する。

> jacac calcPi.java  (これで calcPi.class が生成する)
> java calcPi 1000000  (回数は適宜増減する)

C で書いたプログラム

Ruby はスクリプト言語として高い機能を持っているが，その反面実行速度は犠牲になっている。したがって，モンテカルロ法による積分計算のよう に多数回のループを回す目的には適していない（そもそも高速を要する数値計算では Ruby は全く不利である）。速度の点で優位に立つ言語は C と Fortran である。ここでは，C のソースを紹介しておこう。Cygwin でも gcc を 利用できるので，試してみるとよい。

下のソースは，メルセンヌツイスタを乱数生成に使ったもので，質のよい擬似乱数を高速に作り出して，モンテカルロ法を実行している． コンパイルするにあたっては，MT.h というヘッダファイルを必要とすることに注意．

// calcPi.c
/*
  モンテカルロ法を使って円周率を求める．
  乱数生成には優れた擬似乱数を高速に発生するメルセンヌツイスタのアルゴリズ
  ムを用いる．詳しい情報は同梱の MT.h のコメントを見よ．

コンパイルの仕方：

  MT.hというヘッダファイルをカレントディレクトリにダウンロードしておく．

  次のようにしてコンパイルする

  gcc calcPi.c -lm -O3  -o calcPi

  オプションの意味：
  -lm : 数学ライブラリをリンク
  -O3 : できる限りの最適化を行って高速にする
  -o  : 生成する実行ファイルの名前を指定

  実行するときには次のようにコマンドライン引数を与える

  ./calcPi 1000

*/
#include 
#include 
#include 
#include "MT.h"
#include 
int main(int argc, char *argv[]){
  long i, n_trial, hit = 0;
  double x,y;
  clock_t  time1, time2;
  init_genrand((unsigned)(clock() * UINT_MAX) % ULONG_MAX); // 時刻データで種を蒔く
  n_trial = atol(argv[1]);
  time1 = clock();
  printf("time2= %f\n",(double)time1);
  for(i=0;i < n_trial;i++){
    x = genrand_real1();
    y = genrand_real1();
    if(y < sqrt(1.0-x*x))
      hit ++;
  }
  time2 = clock();
  printf("%12.9f\n",4.0*hit/(double)n_trial);
  printf("経過時間: %7.3f 秒\n", (double)(time2 - time1) / CLOCKS_PER_SEC);
  return 0; /*なくともよい */
}

  gcc calcPi.c -lm -O3  -o calcPi

  calcPi 10000

試してみると，Ruby に比べて圧倒的な速度を体感できるだろう。

モンテカルロ法の誤差

容易に想像できるように，モンテカルロ法による積分計算の計算精度は， かかる時間に比べてよくない。しかし計算が単純で，性質の悪い関数に対しても 安定した結果が得られるという利点はある。 それでは，モンテカルロ法ではどの程度信頼できる値を得られるのだろうか。 それについての考察を次に進めることにする。