無作為抽出を用いることで所望の計算を行う方法を モンテカルロ法 と呼ぶ。理論的な計算から値を求めることが困難な場合などに力を発揮する。最近では、コンピュータ囲碁のプログラムでモンテカルロ法が使用され、棋力が一段と向上した。
(-1)〜(+1)の範囲の乱数を用いて、(x)座標、(y)座標ともにランダムに決めた1000個の点を作る。これを散布図としてプロットすると、以下のようになる。
x <- runif(N, min=-1, max=1)
y <- runif(N, min=-1, max=1)
# アスペクト比1:1で、x軸y軸ともに-1〜+1の範囲をプロット
plot(x, y, asp=1, xlim=c(-1,1), ylim=c(-1,1))
# グリッド線を加える
grid()
ここで、原点からの距離が1未満の点と、距離が1以上になる点とを分類してみよう。
円に入っている点を数えると776個、入っていない点は224個であった。点が散布された正方形の面積は4、半径は(r=1)なので円の面積は(\pi r^2=\pi)である。これは、正方形内部の点（1000個）と円内部の点（776個）の個数の比と等しくなるはずであるから、(\pi\approx 4\times\frac{776}{1000}=3.104)である。
つまり、円周率の近似値は以下のようにして求めることができる。
N <- 500
x <- runif(N, min=-1, max=1)
y <- runif(N, min=-1, max=1)
count <- sum(x*x + y*y < 1)
4 * count / N
## [1] 3.152
上で紹介した、円周率の計算を複数回行ってみよう。以下のプログラムでは一回の計算においてN個の点を用いて円周率を計算し、それを(K)回繰り返している。それぞれの試行の結果をpi.estに貯めておき、最終的にはその平均値とヒストグラムを表示している。
なお、上記の計算とは異なり、第1象限の¼円のみを用いている。
K <- 1000
N <- 100000
pi.est <- rep(0, times=K)
for (k in seq(1,K)) {
x <- runif(N, min=0, max=1)
y <- runif(N, min=0, max=1)
count <- sum(x*x + y*y < 1)
pi.est[k] <- 4*(count / N)
}
cat(sprintf("K=%d N=%d ==> pi=%f\n", K, N, mean(pi.est)))
## K=1000 N=100000 ==> pi=3.141610
hist(pi.est, breaks=50)
rug(pi.est)
中心極限定理により、結果が正規分布に従っている。
あるクイズゲームの優勝者に提示される最終問題。3つのドアがあり、うち1つの後ろには宝が、残り2つにはゴミが置いてあるとする。優勝者は3つのドアから1つを選択するが、そのドアを開ける前にクイズゲームの司会者が残り2つのドアのうち1つを開け、扉の後ろのゴミを見せてくれる。ここで優勝者は自分がすでに選んだドアか、それとも残っているもう1つのドアを改めて選ぶことができる。
さて、ドアの選択を変更することは宝が得られる確率にどの程度影響があるのだろうか。
N <- 10000
prize.at <- floor(runif(N) * 3) + 1 # 宝があるドア (1, 2, or 3)
initial.guess <- floor(runif(N) * 3) + 1 # 最初の選択 (1, 2, or 3)
switch.door <- floor(runif(N) * 2) # ドアを変えるか (1:yes or 0:no)
# ドアを変更して宝が手に入る場合の数を計算
switch.win <- (initial.guess != prize.at) & (switch.door==1)
# ドアを変更せずに宝が手に入る場合の数を計算
stay.win <- (initial.guess == prize.at) & (switch.door==0)
# それぞれの確率を求める
sum(switch.win) / sum(switch.door==1)
## [1] 0.6573777
sum(stay.win) / sum(switch.door==0)
## [1] 0.3369912
確率は約2倍ちがう。つまり、いちど手にしたものは放したくなくなるという「保有バイアス」にあらがって扉の選択を変えることで、2倍の確率で宝を得ることができる。
2の平方根を求める。(x)を0〜2の範囲の一様乱数とし、その2乗（(x)を一辺とする正方形の面積）が2を超えるかどうかを計算する。
N <- 10000
x <- 2 * runif(N)
sum(x^2 < 2) / N * 2
## [1] 1.424
runif()は([0,1))の一様乱数であるため、(x)は(\left[0, 2\right))の範囲となる。すなわち、(x)の値は以下のような性質を持つ。
確率(\sqrt{2}/2)は「(x^2)が2以下の回数」÷「全試行回数」で近似できるので、プログラム中ではsum(x^2 < 2) / N * 2を計算した。