【初心者向け】「モンテカルロ法」概説

[モンテカルロ法  【 Monte Carlo method 】 ](http://e-words.jp/w/%E3%83%A2%E3%83%B3%E3%83%86%E3%82%AB%E3%83%AB%E3%83%AD%E6%B3%95.html)

>ある事象をモデル化した数式や関数があるとき、その定義域に含まれる値をランダムにたくさん生成して実際に計算を行い、得られた結果を統計的に処理することで推定値を得ることができる。数式を解析的に解くのが困難あるいは不可能な場合でも数値的に近似解を求めることができる。
>
>確率論の大数の法則により、試行の回数を増やせば増やすほど解の精度は高まり、無限回の試行を行うと誤差は0に収束することが知られる。
>
>コンピュータでモンテカルロ法の計算を行う場合、値が毎回異なり分布が完全にランダムな真の乱数列を得るには専用のハードウェアが必要になるため、分布の乱雑さは乱数とほぼ変わらないが、一定の計算手順によって確定的に与えられる疑似乱数を用いることが多い。

歴史

>第二次大戦中、米ロスアラモス国立研究所で原子爆弾開発計画に従事していた核物理学者の**スタニスワフ・ウラム**（<font color="blue">Stanislaw Ulam</font>）がこの手法の原型を考案し、同僚の数学者**ジョン・フォン・ノイマン**（<font color="blue">John von Neumann</font>）が戦後、自ら開発したばかりのコンピュータ**ENIAC**で実際に計算を行いその有効性を実証した。
>
>この計算手法は軍事機密としてコードネームが与えられることになったが、二人の同僚の物理学者**ニコラス・メトロポリス**（<font color="blue">Nicholas Metropolis</font>）が、ウラムの叔父がモナコのモンテカルロにあるカジノで借金を抱えた話に着想を得て**モンテカルロ法**（<font color="blue">Monte Carlo method</font> / <font color="blue">Monte Carlo simulation</font>）という名称を考案し、これが採用されて今日までこの名で呼ばれることになった。

##円周率の計算(<font color="blue">全周</font>)

![20190923040544.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597772/0a844a2a-3c70-6438-6d2f-b47e41af2913.png)

```R
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()
```

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

```R
count <- sum(x*x + y*y < 1)
count
>[1] 759
4 * count / N
>3.036
```
色分け対応
![20190922203237.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597772/f1a59be9-ca0c-5dc1-8333-db5c75c7bc80.png)

```R
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()

```
アニメーションさせてみる。
![20190923041709.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597772/300dc4cc-97c6-dc12-9003-56d360995a03.gif)

```R
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")
```

##ヒストグラム（<font color="blue">Histgram</font>）とラグプロット(<font color="blue">Rugplot</font>)の描画
[Rによるヒストグラムの描画](https://data-science.gr.jp/implementation/ida_r_histogram.html)

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

[ラグプロット - RjpWiki](https://data-science.gr.jp/implementation/ida_r_histogram.html)

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

![20190923035812.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597772/906e20a4-e509-641b-918c-18f59075df1a.png)
K=1000 N=100000 ==> pi=3.141853

```R
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)
```
アニメーションさせてみる。
![20190923043618.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597772/9aef6e8e-c4e8-8750-5c1b-236092e3c960.gif)

```R
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")
```

##円周率の計算(<font color="blue">第一象限のみ</font>)

![20190922203405.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597772/dba06fa6-feff-8952-bcb8-035194e8b8f2.png)

```R
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()
```
アニメーションさせてみる。
![20190923042510.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597772/626991d6-ab41-11cd-47fa-00618efde127.gif)

```R
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")
```
ヒストグラムとラグプロット
![20190923040352.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597772/fb85dda1-8d32-54e4-0279-8770324344f9.png)

```R
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
```
アニメーションさせてみる。
![20190923044501.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597772/fa04b19b-5298-c1e9-f336-ccf88042c280.gif)

```R
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")
```

##２の平方根(<font color="blue">sqrt(2)</font>)

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

- <font color="green">x<1</font>である確率は<font color="green">1/2</font>
- <font color="green">x<2</font>である確率は<font color="green">2/2</font>
- <font color="green">sqrt(2)>x</font>である確率は<font color="green">sqrt(2)/2</font>

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

コンセプト

```R
N <- 10000
x <- 2 * runif(N)
sum(x^2 < 2) / N * 2
>[1] 1.4182
```
ヒストグラムとラグプロットのアニメーション
![20190923052104.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597772/ecbde42d-7491-1872-d350-69a88ab08495.gif)

```R
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")
```

##モンティ・ホール問題（<font color="green">Monty Hall problem</font>）。

[モンティ・ホール問題（Monty Hall problem）](https://ja.wikipedia.org/wiki/%E3%83%A2%E3%83%B3%E3%83%86%E3%82%A3%E3%83%BB%E3%83%9B%E3%83%BC%E3%83%AB%E5%95%8F%E9%A1%8C)
>あるクイズゲームの優勝者に提示される最終問題。3つのドアがあり、うち1つの後ろには宝が、残り2つにはゴミが置いてあるとする。優勝者は3つのドアから1つを選択するが、そのドアを開ける前にクイズゲームの司会者が残り2つのドアのうち1つを開け、扉の後ろのゴミを見せてくれる。ここで優勝者は自分がすでに選んだドアか、それとも残っているもう1つのドアを改めて選ぶことができる。さて、ドアの選択を変更することは宝が得られる確率にどの程度影響があるのだろうか。
>
>なんと確率が約2倍ちがう。つまり、いちど手にしたものは放したくなくなるという「**保有バイアス**」にあらがって扉の選択を変えることで、2倍の確率で宝を得ることができる。

コンセプト

```R
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

```

オッズのヒストグラムとラグプロット表示
![20190923062449.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597772/2c3042c5-4d35-e75b-92f0-409e3e131d04.png)

```R
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)
```
**【追伸】**こんな手口もあったのか…
![image.gif](https://cdn-ak.f.st-hatena.com/images/fotolife/o/ochimusha01/20201207/20201207232244.gif)

```R
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")
```
とりあえず、以下続報…