モンテカルロ法による円周率推定（R言語実装）

3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170.............

上の数字は円周率と呼ばれるものであるが，円周率π，円周長l，直径dの時，一番簡単な円周率の定義は，

π = l/d

である．d = 1である円の場合，lがそのまま円周率の値となる．

さて，このπを確率的な方法で”推定”してみる．
一辺が1の正方形に内接する半径r=1の円の扇型を考えてみよう．

正方形自体の面積は1であり，扇型の面積は

(扇型の面積) = (πr^2)/4 = π/4

である．さて，これの面積比を取ってみる．

面積比 = (扇型の面積)/(正方形の面積) = (πr^2)/4 = π/4

次に，この正方形のどこかランダムに点をn個打つとする．そのうち，扇型の中に入った個数をmとする．このランダムに点を打つとき，扇型の中に点が打たれる確率は，

(扇型の中に点が打たれる確率) = m/n

である．この確率は，先に出した面積比と等号で結ぶことができる．

(面積比) = (扇型の中に点が打たれる確率)
π/4 = m/n
∴π = 4m/n

つまり，扇型の中に点が打たれる確率で円周率を推定することができる．この方法を「モンテカルロ法」とよび，数値積分などの計算にも利用されている．

さて，ではR言語を用いてモンテカルロ法による円周率の推定を行なってみよう．次のプログラムを用いて推定を行った．

mcPi <- function(points){
  inR <- 0 #扇型に入った数
  cols <- NULL #plot色分け用
  range <- c(0,1) #plot時の描画範囲
  x <- runif(points,0,1) #一様乱数でx座標を決定
  y <- runif(points,0,1) #一様乱数でy座標を決定
  r <- x^2+y^2 #原点から各点がどのくらい離れているか
  for(i in r){
    if(i<=1){ #もし，原点からの距離が1以下であれば扇型の中
      inR <- inR + 1
      cols <- c(cols,"red") #扇型内であれば赤色でプロット
    }else{
      cols <- c(cols,"blue") #扇型外であれば青色でプロット
    }
  }
  return_data <- 4*inR/points #円周率計算
  plot(x,y,col=cols,xlim=range,ylim=range,pch=20,main="Estimate pi value with Monte Carlo method.",sub=sprintf("%s points. pi = %1.6f",points, return_data))
  return(return_data)
}

このmcPi()は引数に応じて点を打っていき，円周率を計算するプログラムである．

実行例：
mcPi(10000) #10000回点を打って円周率を推定

結果例：
[1] 3.1516 ←これが推定された円周率の値．毎回違う値になる．

次のような画像も一緒に出力される．

当然のことだが，打つ点が多ければ多いほどよい精度になりやすい（πの値により近づきやすくなる）．