Montecarlo Method

Montecarlo法を用いて円周率\(\pi\)の近似値を求める。

　

半径が\(r\)の円と、その円に外接する正方形の面積を考える。 円の面積は\(\pi r^{2}\)であり、正方形の面積は\(4r^{2}\)である。

これらの面積比を\(a\)として計算すると、\(a=\dfrac{\pi r^{2}}{4r^{2}}=\dfrac{\pi}{4}\)となる。

つまり、\(a\)の値を推定することで円周率を推定できる。

ここで、この正方形内にランダムに\(n\)個の点を打ち、そのうち円内部に入った点の個数を\(m\)とする。

このとき、円内部に点が打たれる確率は\(\dfrac{m}{n}\)である。

よって、\(\dfrac{\pi}{4}=\dfrac{m}{n}\)すなわち\(\pi=\dfrac{4m}{n}\)となる。

これによって円周率を推定することができる。

以下のプログラムでは簡単のため、円の半径を\(1\)とし、第一象限のみで考えるものとする。

montecarlo <- function(n){
  m <- 0 #円内部にある点の個数; 始めは0個とする
  cols <- NULL #plotする点に色をつけるためのもの
  range <- c(0,1) #範囲を第一象限とする
  x <- runif(n,0,1) #一様乱数でx座標を定める; runif(試行回数,開始値,終了値)
  y <- runif(n,0,1) #一様乱数でy座標を定める
  r <- x^2+y^2 #円の範囲
  for(i in r){ 
    if(i<=1){
      m <- m +1 #円内部であれば、円内部の点の個数を1つ増やす
      cols <- c(cols, "red") #更に赤色でプロットする
    } else{
      cols <- c(cols,"blue") #円外部であれば青色でプロットする
    }
  }
# forループは以下のようにも書ける
#  inner <- (r<=1)
#  m <- sum(inner)
#  cols <- ifelse(inner,"red","blue")
  return_data <- 4*m/n #円周率の計算
  plot(x,y,col=cols,xlim=range,ylim=range,asp=1)
  #点をプロットする,xとyの範囲は上で定めたもの,x軸とy軸の比率を1にする
  return(return_data)
}

\(n=1000\)として結果を出力すると、図と出力結果は以下である。

## [1] 3.088

また、\(n=10000\)として結果を出力すると結果は以下である。

## [1] 3.1456