1 円周率をシミュレーションで計算

1.1 解説

シミュレーションから円周率を求める方法を説明します（有名な例）。求め方は簡単です。まず一辺の長さが1となる正方形を書きます。次に，半径1の円の方程式は，\(x^2 + y^2 =1\)ですので， \(y\)を\(x\)の関数として表すと， \[ y = \sqrt{1-x^2} \] となります。図示すると以下のようになります。

a <- seq(from = 0, to = 1, length=100)
b <- sqrt( 1 - a^2)

plot(a, b, type="l")
lines(c(1,0),c(0,0))
lines(c(0,0),c(0,1))
lines(c(1,0),c(1,1))
lines(c(1,1),c(0,1))

ここでは，seq()を用いて0から1までの値をとる100の数からなる等差数列を作成しxとする。 \(\sqrt{1 - x^2}\)としたものをyとする。

1.2 乱数

n = 100 # サンプルサイズ
x <- runif(n,0,1)
y <- runif(n,0,1)

1.3 シミュレーション

plot(a, b, type="l")
lines(c(1,0),c(0,0))
lines(c(0,0),c(0,1))
lines(c(1,0),c(1,1))
lines(c(1,1),c(0,1))
inner <- 0
for (i in 1:n) { # iに1から100まで以下を繰り返す
        points(x[i], y[i], pch=4, col="red") # 座標 (x_i, y_i)の赤色で×印
            if ( x[i]^2 + y[i]^2 <= 1 ){ # x^2+y^2<=1なら，
            inner <- inner + 1 # innerに1を加えていく。
            }
    }

1.4 結果

ratio <- inner/n*4
ratio

## [1] 3.04

1.5 試行を増やしてみる。

1000個のごまをばらまく。

n = 1000 # サンプルサイズ
x <- runif(n,0,1)
y <- runif(n,0,1)
inner1K <- 0
for (i in 1:n) {
            if ( x[i]^2 + y[i]^2 <= 1 ){
            inner1K <- inner1K + 1
            }
    }

何度も書くのが面倒なので，関数にしておく。

enshu <- function(n = 100) {
  trial = n
  x <- runif(n,0,1)
  y <- runif(n,0,1)
  inner <- 0
  for (i in 1:n) {
            if ( x[i]^2 + y[i]^2 <= 1 ){
            inner <- inner + 1
            }
  }
  ratio <- inner/n*4
  return(ratio)
}

1000個のごまをばらまく。

res1K <- enshu(1000)

10000個のごまをばらまく

res10K  <- enshu(10000)
res100K <- enshu(100000)
res1M <- enshu(1000000)

結果を比べる

print(c(ratio,res10K,res100K,res1M))

## [1] 3.04000 3.12800 3.14188 3.14044

となり，試行回数が増えるほど円周率\(3.141592658589 \cdots\)により近い結果となっていることが分かります。

N = 100 for ( i in 1:N)

円周率の求め方