モンテカルロ法【 Monte Carlo method 】

ある事象をモデル化した数式や関数があるとき、その定義域に含まれる値をランダムにたくさん生成して実際に計算を行い、得られた結果を統計的に処理することで推定値を得ることができる。数式を解析的に解くのが困難あるいは不可能な場合でも数値的に近似解を求めることができる。

確率論の大数の法則により、試行の回数を増やせば増やすほど解の精度は高まり、無限回の試行を行うと誤差は0に収束することが知られる。

コンピュータでモンテカルロ法の計算を行う場合、値が毎回異なり分布が完全にランダムな真の乱数列を得るには専用のハードウェアが必要になるため、分布の乱雑さは乱数とほぼ変わらないが、一定の計算手順によって確定的に与えられる疑似乱数を用いることが多い。

歴史

第二次大戦中、米ロスアラモス国立研究所で原子爆弾開発計画に従事していた核物理学者のスタニスワフ・ウラム（Stanislaw Ulam）がこの手法の原型を考案し、同僚の数学者ジョン・フォン・ノイマン（John von Neumann）が戦後、自ら開発したばかりのコンピュータENIACで実際に計算を行いその有効性を実証した。

この計算手法は軍事機密としてコードネームが与えられることになったが、二人の同僚の物理学者ニコラス・メトロポリス（Nicholas Metropolis）が、ウラムの叔父がモナコのモンテカルロにあるカジノで借金を抱えた話に着想を得てモンテカルロ法（Monte Carlo method / Monte Carlo simulation）という名称を考案し、これが採用されて今日までこの名で呼ばれることになった。

##円周率の計算(全周)

N<-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), main="Monte Carlo simulation",xlab="Real Expanse", ylab="Imaginal Expanse")

par(new=T)#上書き指定
theta00<- seq(0, -2*pi, length=360)
theta_cos<-cos(theta00)
theta_sin<-sin(theta00)
plot(theta_cos,theta_sin,asp=1,xlim=c(-1,1), ylim=c(-1,1), type="l", main="",xlab="", ylab="")

# グリッド線を加える
grid()

点が散布された正方形の面積は2*2=4
半径r=1なので円の面積はπ*1^2=π
πの値は正方形内部の点（1000個）と円内部の点の個数（count）の比と等しくなるはずである。従ってπ≈4*count/Nとなる。

count <- sum(x*x + y*y < 1)
count
>[1] 759
4 * count / N
>3.036

色分け対応

N<-1000
x0 <- runif(N, min=-1, max=1)
y0 <- runif(N, min=-1, max=1)
x_in<-c(NULL)
y_in<-c(NULL)
x_out<-c(NULL)
y_out<-c(NULL)
for (i in 1:N){
if(x0[i]*x0[i]+y0[i]*y0[i]<1)
{x_in<-c(x0[i],x_in);y_in<-c(y0[i],y_in)}
else
{x_out<-c(x0[i],x_out);y_out<-c(y0[i],y_out)}
}

# アスペクト比1:1で、x軸y軸ともに-1〜+1の範囲をプロット
plot(x_in, y_in, asp=1,col=rgb(1,0,0), xlim=c(-1,1), ylim=c(-1,1), main="Monte Carlo simulation",xlab="Real Expanse", ylab="Imaginal Expanse")

par(new=T)#上書き指定
plot(x_out,y_out,asp=1,col=rgb(0,0,1),xlim=c(-1,1), ylim=c(-1,1),main="",xlab="", ylab="")

# グリッド線を加える
grid()

アニメーションさせてみる。

Monte_Carlo_simulation01<-function(x){
#xの値(試行回数)は実際には使わない。

N<-1000
x0 <- runif(N, min=-1, max=1)
y0 <- runif(N, min=-1, max=1)
x_in<-c(NULL)
y_in<-c(NULL)
x_out<-c(NULL)
y_out<-c(NULL)
for (i in 1:N){
if(x0[i]*x0[i]+y0[i]*y0[i]<1)
{x_in<-c(x0[i],x_in);y_in<-c(y0[i],y_in)}
else
{x_out<-c(x0[i],x_out);y_out<-c(y0[i],y_out)}
}

# アスペクト比1:1で、x軸y軸ともに-1〜+1の範囲をプロット
plot(x_in, y_in, asp=1,col=rgb(1,0,0), xlim=c(-1,1), ylim=c(-1,1), main="Monte Carlo simulation",xlab="Real Expanse", ylab="Imaginal Expanse")

par(new=T)#上書き指定
plot(x_out,y_out,asp=1,col=rgb(0,0,1),xlim=c(-1,1), ylim=c(-1,1),main="",xlab="", ylab="")

# グリッド線を加える
grid()

}

library("animation")
Time_Code=c(1,15,30,45,60,75,90,105,120,135,150,165,180,195,210,225,240,255,270,285,300,315,330,345)
saveGIF({
for (i in 1:16){
　　Monte_Carlo_simulation01(i)
}
}, interval = 0.1, movie.name = "TEST001.gif")

##ヒストグラム（Histgram）とラグプロット(Rugplot)の描画
Rによるヒストグラムの描画

データの階級数を客観的に決定するための指標として，スタージェスの公式があり、データの個数nと階級数kをk≒1+log(N,base=2)のように決定する．コマンドhist はデフォルトではこの式に基づいて階級数および階級幅を決定している．

ラグプロット - RjpWiki

ラグ(敷物)プロットは一次元データをバーコード風に表現する。単独で使うものではないが、他のプロットの x,y 軸の飾りとして使うとなかなか見栄えがよい。

K=1000 N=100000 ==> pi=3.141853

K <- 1000
N <- 100000
pi.est <- c(NULL)
for (k in seq(1,K)) {
x <- runif(N, min=-1, max=1)
y <- runif(N, min=-1, max=1)
count <- sum(x*x + y*y < 1)
pi.est<-c(4*(count/N),pi.est)
}
cat(sprintf("K=%d N=%d ==> pi=%f\n", K, N, mean(pi.est)))
## K=1000 N=100000 ==> pi=3.141659
hist(pi.est, breaks=50)
rug(pi.est)

アニメーションさせてみる。

MC_Hist01<-function(x){
#xの値(試行回数)は実際には使わない。

K <- 1000
N <- 100000
pi.est <- c(NULL)
for (k in seq(1,K)) {
x <- runif(N, min=-1, max=1)
y <- runif(N, min=-1, max=1)
count <- sum(x*x + y*y < 1)
pi.est<-c(4*(count/N),pi.est)
}
hist(pi.est, breaks=50)
rug(pi.est)

}

library("animation")
saveGIF({
for (i in 1:16){
　　MC_Hist01(i)
}
}, interval = 0.1, movie.name = "MCH001.gif")

##円周率の計算(第一象限のみ)

N<-1000
x0 <- runif(N, min=0, max=1)
y0 <- runif(N, min=-0, max=1)
x_in<-c(NULL)
y_in<-c(NULL)
x_out<-c(NULL)
y_out<-c(NULL)
for (i in 1:N){
if(x0[i]*x0[i]+y0[i]*y0[i]<1)
{x_in<-c(x0[i],x_in);y_in<-c(y0[i],y_in)}
else
{x_out<-c(x0[i],x_out);y_out<-c(y0[i],y_out)}
}

# アスペクト比1:1で、x軸y軸ともに-1〜+1の範囲をプロット
plot(x_in, y_in, asp=1,col=rgb(1,0,0), xlim=c(0,1), ylim=c(0,1), main="Monte Carlo simulation",xlab="Real Expanse", ylab="Imaginal Expanse")

par(new=T)#上書き指定
plot(x_out,y_out,asp=1,col=rgb(0,0,1),xlim=c(0,1), ylim=c(0,1),main="",xlab="", ylab="")

# グリッド線を加える
grid()

アニメーションさせてみる。

Monte_Carlo_simulation02<-function(x){
#xの値(試行回数)は実際には使わない。

N<-1000
x0 <- runif(N, min=0, max=1)
y0 <- runif(N, min=0, max=1)
x_in<-c(NULL)
y_in<-c(NULL)
x_out<-c(NULL)
y_out<-c(NULL)
for (i in 1:N){
if(x0[i]*x0[i]+y0[i]*y0[i]<1)
{x_in<-c(x0[i],x_in);y_in<-c(y0[i],y_in)}
else
{x_out<-c(x0[i],x_out);y_out<-c(y0[i],y_out)}
}

# アスペクト比1:1で、x軸y軸ともに-1〜+1の範囲をプロット
plot(x_in, y_in, asp=1,col=rgb(1,0,0), xlim=c(0,1), ylim=c(0,1), main="Monte Carlo simulation",xlab="Real Expanse", ylab="Imaginal Expanse")

par(new=T)#上書き指定
plot(x_out,y_out,asp=1,col=rgb(0,0,1),xlim=c(0,1), ylim=c(0,1),main="",xlab="", ylab="")

# グリッド線を加える
grid()

}

library("animation")
Time_Code=c(1,15,30,45,60,75,90,105,120,135,150,165,180,195,210,225,240,255,270,285,300,315,330,345)
saveGIF({
for (i in 1:16){
　　Monte_Carlo_simulation02(i)
}
}, interval = 0.1, movie.name = "TEST002.gif")

ヒストグラムとラグプロット

K <- 1000
N <- 100000
pi.est <- c(NULL)
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<-c(4*(count/N),pi.est)
}
cat(sprintf("K=%d N=%d ==> pi=%f\n", K, N, mean(pi.est)))
## K=1000 N=100000 ==> pi=3.141659
hist(pi.est, breaks=50)
rug(pi.est)

結果例
K=1000 N=100000 ==> pi=3.141352

アニメーションさせてみる。

MC_Hist02<-function(x){
#xの値(試行回数)は実際には使わない。

K <- 1000
N <- 100000
pi.est <- c(NULL)
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<-c(4*(count/N),pi.est)
}
hist(pi.est, breaks=50)
rug(pi.est)

}

library("animation")
saveGIF({
for (i in 1:16){
　　MC_Hist02(i)
}
}, interval = 0.1, movie.name = "MCH002.gif")

##２の平方根(sqrt(2))

xを範囲0〜2の一様乱数とし、その2乗（を一辺とする正方形の面積x）が2を超えるかどうかを計算する。runif()は[0,1)の一様乱数であるから、x*2は[0,2)の範囲となる。すなわち、xの値は以下のような性質を持つ。

x<1である確率は1/2
x<2である確率は2/2
sqrt(2)>xである確率はsqrt(2)/2

確率sqrt(2)/2は「x^2が2以下の回数÷全試行回数」で近似できるので、プログラム中ではsum(x^2<2)/N*2を計算。

コンセプト

N <- 10000
x <- 2 * runif(N)
sum(x^2 < 2) / N * 2
>[1] 1.4182

ヒストグラムとラグプロットのアニメーション

Sqrt2_01<-function(x){
#xの値(試行回数)は実際には使わない。
K <- 1000
N <- 100000
sqrt2.est <- c(NULL)
for (k in seq(1,K)) {

N <- 10000
x <- 2 * runif(N)
sqrt2.est<-c(sum(x^2 < 2) / N * 2,sqrt2.est)
}
hist(sqrt2.est, breaks=50)
rug(sqrt2.est)
}
library("animation")
saveGIF({
for (i in 1:16){
　　Sqrt2_01(i)
}
}, interval = 0.1, movie.name = "Sqrt2_01.gif")

##モンティ・ホール問題（Monty Hall problem）。

モンティ・ホール問題（Monty Hall problem）

あるクイズゲームの優勝者に提示される最終問題。3つのドアがあり、うち1つの後ろには宝が、残り2つにはゴミが置いてあるとする。優勝者は3つのドアから1つを選択するが、そのドアを開ける前にクイズゲームの司会者が残り2つのドアのうち1つを開け、扉の後ろのゴミを見せてくれる。ここで優勝者は自分がすでに選んだドアか、それとも残っているもう1つのドアを改めて選ぶことができる。さて、ドアの選択を変更することは宝が得られる確率にどの程度影響があるのだろうか。

なんと確率が約2倍ちがう。つまり、いちど手にしたものは放したくなくなるという「保有バイアス」にあらがって扉の選択を変えることで、2倍の確率で宝を得ることができる。

コンセプト

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.664993
sum(stay.win) / sum(switch.door==0)
## [1] 0.3393746

結果例（オッズとしては前者が後者の２倍）

[1] 0.6675835
[1] 0.3266802

オッズのヒストグラムとラグプロット表示

MHP.est<-c(NULL)
for(i in 1:1000){
#iの値(試行回数)は実際には使わない。
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)
# それぞれの確率を求める
a<-sum(switch.win) / sum(switch.door==1)
b<-sum(stay.win) / sum(switch.door==0)
#オッズを記録
MHP.est<-c(a/b,MHP.est)

}
hist(MHP.est, breaks=50)
rug(MHP.est.est)

**【追伸】**こんな手口もあったのか…

circle_check<-function(x,y){
ifelse(sqrt(x^2+y^2)<=1,TRUE,FALSE)
}

Monte_Carlo_simulation02<-function(x){
#xの値(試行回数)は実際には使わない。

N<-1000
x0 <- runif(N, min=-1, max=1)
y0 <- runif(N, min=-1, max=1)

# アスペクト比1:1で、x軸y軸ともに-1〜+1の範囲をプロット
plot(x0,y0,asp=1,xlim=c(-1,1),ylim=c(-1,1), main="Monte Carlo simulation",xlab="Real Expanse", ylab="Imaginal Expanse",col=ifelse(circle_check(x0,y0),rgb(1,0,0),rgb(0,0,1)))

# グリッド線を加える
grid()
}

library("animation")
Time_Code=c(1,15,30,45,60,75,90,105,120,135,150,165,180,195,210,225,240,255,270,285,300,315,330,345)
saveGIF({
for (i in 1:16){
　　Monte_Carlo_simulation02(i)
}
}, interval = 0.1, movie.name = "TEST001.gif")

とりあえず、以下続報…