モンテカルロ法の誤差を考える

概要

PIcalc

円周率をモンテカルロ法で求める問題では， 図の正方形の部分に一様乱数をランダムに打ち込んだ場合， 円弧の中に入る確率は$p = \pi/4 (= 0.7853981633974483\ldots) $になる。 打ち込んだ乱数の数を $n$, 入った数を $x$ とすると， 円周率　$\pi$の近似値として $4x/n$ をとればよいことは，ちょっと計算してみればわかるだろう。

さて， 十分に多数の乱数を使えば精度よく円周率が求まるだろうというのは自然な考えだ。 それでは，どれくらい多数にすれば，どれくらい精度が高まるのだろうか？ それを知っておかないと，結果に自信をもつことはできないよね。

確率論から $x$の分布は二項分布になることがわかり，さらに $x$ の 期待値は $\mu = np$，標準偏差は $\sigma = \sqrt{np(1-p)}$ である。 この標準偏差から，精度を求めることができる。

標準偏差は $x$ の広がりの幅を決めるパラメータだ。$n$ が大きければ $x$ の分布は 正規分布とみなせるので，正規分布表や 正規分布のパーセント点を参照すれば，99% の乱数が円弧に入る範囲を見積もることができる． 実際，$x$が平均値を中心にして $\pm 2.576 \sigma$ の範囲に入ることが，表から分かる。 上に述べたように， $x$ に $4/n$を掛ければ $\pi$の推定値がえられるから， $4x/n$ の値の 99%が入る範囲の幅は次のようになる。

\[ \pm 2.576 \sigma \times \frac4n = \frac{\pm 4\times2.576\sqrt{np(1-p)} }{n} = \frac{\pm4.230}{\sqrt{n}} \]

一方，分布の中央は \[ \mu \times \frac4n = 4p \] であるから，結局 $4x/n$の値は次の区間内に 99% の確率で集まることになる。 \[ \left[ 4p - \frac{4.230}{\sqrt{n}}, 4p + \frac{4.230}{\sqrt{n}} \right] \]

実験

それでは，モンテカルロ法で得られる $\pi$ の値の分布を調べるための実験を行っていくことにしよう。高速で実行できるCのプログラムと，テキストやファイルの扱いに優れた Ruby のスクリプト を併用すると効率がよい。

C のプログラム

Cのプログラムの名前は calcPi2.c，それをコンパイルした実行プログラムは calcPi2 とする。ただし Windows の場合，バイナリの実行プログラムの拡張子は .exe とすることになっているので， calcPi2.exeになる。このプログラムは， コマンドライン引数で 「打ち込む」乱数の個数（統計学的には「標本の大きさ」である）（ここでは 10000）を与えて実行する。すると瞬時に$\pi$の推定値の計算結果 が標準出力に返される。 

> ./calcPi2 10000 
   3.157600000
これをたとえば 1000回やって，たくさんの推定値を集める。その分布と，標本の大きさとの関係を調べ，理論値と比較してみようというのが，この実験の流れである。

一回の計算で「打ち込む」乱数の数を n_sampleとしよう。

プログラムを走らせてみる。

計算のための C プログラム

次のプログラムは，すでに紹介したプログラムから，処理時間の出力を省いて，$\pi$ の推定値だけを出力するようにしたものである。 

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

 コンパイルは次のようにする
   gcc calcPi2.c -lm -O3  -o calcPi2

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

   実行するときには次のようにコマンドライン引数を与える
   ./calcPi2 1000
*/
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <limits.h>
#include "MT.h"
int main(int argc, char *argv[]){
  long i, n_trial, hit = 0;
  double x,y;
  // 精密な時刻データを使って種を蒔く
  init_genrand((unsigned)(clock() * UINT_MAX) % ULONG_MAX); 
  n_trial = atol(argv[1]);
  for(i=0;i < n_trial;i++){
    x = genrand_real1();
    y = genrand_real1();
    if(y < sqrt(1.0-x*x))
      hit ++;
  }
  printf("%12.9f\n",4.0*hit/(double)n_trial);
  return 0; 
}