library(tidyverse)
library(ggforce)

1. モンテカルロ・シミュレーションとは

モンテカルロ法 (Monte Carlo method)

  • 無作為に抽出された乱数を用い、数値計算やシミュレーションを行う手法
  • モンテカルロ法を用いたシミュレーションがモンテカルロ・シミュレーション
  • モンテカルロはヨーロッパのモナコ公国内の一つの地区であり、カジノで有名

  • カジノではサイコロやルーレット、無作為に配られたカードや乱数が頻繁に使われることが名前の由来
  • 計算があまりにも複雑だったり、実質的に代数で解が得られない解析学上の問題などに強みを持つ手法
  • パソコンの性能が著しく向上
    →　モンテカルロ法（マルコフ連鎖モンテカルロ法; MCMC）によって、ようやくベイズ統計学が使えるようになった
  • モンテカルロ法から得られた結果には常に誤差が存在する
  • 特に生成された乱数が少ない（= 試行回数が少ない）場合
    →　この誤差は大きくなる
  • 近年のパソコンの性能が飛躍的に発達
    →　かなり小さな誤差で正確な結果が得られるようになった

1.1 乱数の生成

sample() 関数によるサンプリング

  • 無作為に値を抽出する 2 つの方法：
  1. 値の集合から無作為に値を抽出する方法　　
  2. 正規分布などの確率分布から値を抽出する方法
(1) 値の集合から無作為に値を抽出する方法
  • sample()関数を用い、値の集合から無作為に値を抽出する方法

sample() 関数によるサンプリング sample(x = 値の集合ベクトル,
　　　size = 抽出の回数,
　　　replace = 復元抽出の有無,
　　　prob = 各要素が抽出される確率)

  • 例えば、x は 1 から 10 までの値を含む集合ベクトルとする　　
x <- c(1, 2, 3, 4, 5, 6)
  • ｘから無作為に 5 つの値を選びたければ
sample(x,
       size = 5)
[1] 6 2 4 5 3
  • この場合、一度抽出された要素は、二度と抽出されない
    →　何度やっても同じ値が選ばれることはない
  • しかし、復元抽出の有無を指定する引数 (replace)を replace = TRUE と指定
    →　選ばれた数字を x に戻して、数字を選ぶ
    →　同じ値が選ばれる可能性がでてくる
sample(x,
       size = 5, 
       replace = TRUE)
[1] 1 6 3 1 3
  • prob は各要素が抽出される確率
  • x の実引数と同じ長さの numeric 型ベクトルを指定
  • prob の実引数の総和は 1 であることが望ましい
  • 総和が1でない場合　→　自動的に総和が 1 になるよう正則化を行う
  • 例えば、4, 5, 6 が 1, 2, 3 より 2 倍の確率で選ばれるように指定したければ、次のように指定する
sample(x,
       size = 5, 
       replace = TRUE,
       prob = c(1, 1, 1, 2, 2, 2))
[1] 4 5 4 4 4

  • つまり次の二つのコマンドは同じ意味

  • prob = c(1, 1, 1, 2, 2, 2)

  • prob = c(1/9, 1/9, 1/, 2/9, 2/9, 2/9)

  • サイコロを 3 回振るコードを書くなら、抽出の回数 size = 3

  • サイコロであれば、一回出た目も抽出される可能性がある
    →　復元抽出を行う必要あり（replace = TRUE）

  • 各目が出る確率は1/6

  • 各値が抽出される確率が等しいので、省略可能

sample(1:6,
       3,
       replace = TRUE)
[1] 2 4 3
  • サイコロに仕掛けをして、奇数の目 (1, 3, 5) が偶数の目 (2, 4, 6)よりも2倍の確率で出るように細工をしたとすれば
sample(x,
       size = 5, 
       replace = TRUE,
       prob = c(2, 1, 2, 1, 2, 1))
[1] 1 4 1 5 2
  • このサイコロを1万回振り、それぞれの目が出た回数を棒グラフとして示す
  • 無作為に抽出された値であれば、各目が出る回数は等しいはず
  • table() 関数を使うと、各要素が出現した回数が出力される
    →　このオブジェクトを barplot() 関数に渡すと棒グラフを作成できる
dice <- sample(1:6, 
               10^4, 
               replace = TRUE,
               prob = c(2, 1, 2, 1, 2, 1))
table(dice) |> 
  barplot()

  • サイコロを1万回振ったシミュレーションでも、奇数の目 (1, 3, 5) が偶数の目 (2, 4, 6) よりも2倍の確率で出ることがわかる

1.2 確率分布からの乱数生成（主要な分布）

関数名 確率分布 パラメータ
rbeta() ベータ分布 n, shape1, shape2
rbinom() 二項分布 n, size, prob
rcauchy() コーシー分布 n, location, scale
rchisq() \(X^2\)分布 n, df
rexp() 指数分布 n, rate
rf() F 分布 n,
rgamma() ガンマ分布 n, shape, scale
rgeom() 幾何分布 n, prob
rhyper() 超幾何分布 nn, m, n, k
rlnorm() 対数正規分布 n, meanlog, sdlog
rmultinom() 多項分布 n, size, prob
rnbinom() 負の二項分布 n, size, prob
rnorm() 正規分布 n, mean, sd
rpois() ポアソン分布 n, lambda
rt() t分布 n, df
runif() 一様分布 n, min, max
rweibull() ワイブル分布 n, shape, scale
mvtnorm::rmvnorm() 多変量正規分布 n, mean, sigma

1.3 シードについて

  • 特定の分布から乱数を抽出する場合、抽出の度に値が変わる
  • 例えば、平均0、標準偏差1の正規分布（標準正規分布）から 3 つの値を抽出し、小数点 2 桁に丸める作業を 3 回繰り返す
rnorm(3) |> 
        round(2)
[1] 0.49 0.76 0.84
rnorm(3) |> 
        round(2)
[1] -0.64 -0.51 -0.88
rnorm(3) |> 
        round(2)
[1] -2.37 -0.12 -1.98
  • 抽出の度に結果が変わる
  • 1 回きりのシミュレーションではこれで問題ない
  • しかし、同じシミュレーションから同じ結果を得るためには、乱数を固定する必要がある
  • このような時に使うのがシード（seed）
  • シードが同じなら抽出される乱数は同じ値
  • シードの指定は set.seed(numeric型スカラー)
  • 例えば、シードを 20220826 にして実行してみる
set.seed(20220826)
rnorm(3) |> 
        round(2)
[1] -0.14 -0.87 -0.29
set.seed(20220826)
rnorm(3) |> 
        round(2)
[1] -0.14 -0.87 -0.29
  • シードを指定すると乱数が固定され、シミュレーション結果の再現ができる

2. Birthday Paradox

  • 誕生日問題 (Birthday Paradox)
  • 「ある集団の中で、誕生日が同じペアが一組以上いるためには、何人必要か？」という問題
  • 2月29日を除外すると、365人いれば 100% の確率で同じ誕生日の人が一組いる
  • では、この人数が半分（365/2）なら、同じ誕生日のペアが一組いる確率は 50% になるのか？

正解 23 人いれば、同じ誕生日のペアが一組いる確率が 50% になる

R を使ったシミュレーションを使って確認してみる

  • 学生の数を n_student と設定する
  • 1 から 365 までの公差 1 の等差数列から、n_student個の値を復元抽出する
  • ここでは 10 人の学生がいるとする
n_student <- 10
Birth_vec <- sample(1:365,
                    n_student,
                    replace = TRUE)
Birth_vec
 [1] 140  55 363 170 292 335 345   6 179   6
  • ここで数字が 1 組重複すれば　→　誕生日が同じ学生が 1 組いる
  • 重複する要素（＝数字）があるかどうかを確認するには duplicated()関数を使う
  • ある要素が他の要素と重複するなら TRUE が、なければ FALSE が表示される
duplicated(Birth_vec)
 [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE
  • 1 つでも TRUE が含まれていれば　→　同じ誕生日の人が 2 人以上はいるという意味
  • 1 つでも TRUE が含まれているかどうかを確かめるにはany()関数を使う
    →　1 つでも TRUEがあれば（= 同じ誕生日のペアがいれば）→　返り値は TRUE
    →　全て FALSE なら　→　返り値は FALSE
any(duplicated(Birth_vec))
[1] TRUE
  • ここでは 6 が重複している　→　 TRUE と表示
  • この作業を100回繰り返してみる
n_student <- 10   # 学生数
n_trials  <- 100  # 試行回数


Result_vec <- rep(NA, n_trials) # 結果を格納する空ベクトル Result_vec を用意

# 反復処理
for (i in 1:n_trials) {
    Birth_vec     <- sample(1:365, n_student, replace = TRUE)
    Result_vec[i] <- any(duplicated(Birth_vec))
}

Result_vec
  [1] FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE
 [13] FALSE  TRUE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE
 [25] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [37] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [49] FALSE  TRUE FALSE  TRUE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE  TRUE
 [61] FALSE FALSE  TRUE FALSE  TRUE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE
 [73] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [85] FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE
 [97] FALSE FALSE FALSE FALSE
sum(Result_vec)
[1] 12
  • TRUE が 12 個　→　100 人のうち 12 人（12%）の誕生日が同じ　　
  • この繰り返しを 1 万回繰り返す　→　確率に近似できそう
n_student <- 10    # 学生数
n_trials  <- 1000 # 試行回数

# 結果を格納する空ベクトルを用意する
Result_vec <- rep(NA, n_trials)

# 反復処理
for (i in 1:n_trials) {
    Birth_vec     <- sample(1:365, n_student, replace = TRUE)
    Result_vec[i] <- any(duplicated(Birth_vec))
}

sum(Result_vec)
[1] 110
  • 1万回の試行から TRUE が出た回数は 1190回（= 11.9%）
  • つまり、11人集まれば誕生日が同じペアが 2 つ以上の確率は約 11.9%ということ

2.1 シミュレーションと理論値の比較

2.1.1 シミュレーション（試行回数 = 100）

  • ここでは試行回数を 100 に固定
  • 学生数を 2 から 100 まで調整しながら、同じ誕生日の人が 2 人以上いる割合をシミュレーションしてみる
Prob_df <- tibble(Students = 2:100, # 学生数を 2 から 100 まで調整  
                  Probs    = NA) # Probs は空に指定

for (i in 1:nrow(Prob_df)) {
    
    Result_vec <- rep(NA, 100) # 試行回数を 100 に固定 
    
    for (j in 1:100) {        # 試行回数を 100 に固定 
        Birth_vec     <- sample(1:365, Prob_df$Students[i], replace = TRUE)
        Result_vec[j] <- any(duplicated(Birth_vec))
    }
    
    Prob_df$Probs[i] <- mean(Result_vec)
}
  • 結果を確認する
DT::datatable(Prob_df)
  • 生徒数が 2 の時に、この 2 人の誕生日が同じ確率は 0
  • 同じ誕生日の人が 2 人以上いる割合 (Probs) が 0.5 (= 50%) になる人数 (Students) を確かめる
Prob_df |> 
  filter(Probs >= 0.4 & Probs <= 0.6)
# A tibble: 7 × 2
  Students Probs
     <int> <dbl>
1       19  0.43
2       20  0.43
3       21  0.5 
4       22  0.51
5       23  0.59
6       24  0.58
7       25  0.59
  • 22人か23人のあたりだとわかる
  • x 軸を「学生数」、y 軸を「同じ誕生日の人が2人以上いる割合 (%)」と指定して結果を可視化みる
Prob_df %>%
    ggplot() +
    geom_line(aes(x = Students, 
                  y = Probs * 100), 
              size = 1) +
    geom_hline(yintercept = 50, # 50%に赤い点線を引く  
               color = "red", 
               linetype = 2) +
    labs(x = "学生数",
         y = "同じ誕生日の人が2人以上いる割合 (%)\n(試行回数 = 100)") +
    theme_minimal(base_family = "HiraKakuProN-W3")

シミュレーション結果（試行回数 = 100） 50人ほどいれば、ほぼ確実に同じ誕生日のペアが一組以上いる

2.1.2 理論値

  • n 人いる場合、同じ誕生日の人が 2 人以上いる確率 \(P(n)\)

\[P(n) = 1 - \frac{365!}{365^n(365-n)!}\]

  • \(3! = 3 × 2 × 1 = 6\)
    →　\(365! = 365 × 364 ×....× 1\)

  • R では factorial() 関数を使って計算できる

  • 50 人いるときに同じ誕生日のペアがいる確率は

1 - factorial(365) / (365^50 * factorial(365 - 50))
[1] NaN
  • NaN は「計算不可」という意味
  • 数字が大きすぎて計算できない

NaNについて - \(100!\)を計算してみる

factorial(100)
[1] 9.332622e+157

  • 小数点以下に 0 が 157 個もつくほど大きな値

  • \(200!\)を計算してみる

factorial(200)
[1] Inf
  • 値が大きすぎて計算不可
  • 上の式を計算可能な式に変形する

\[P(n) = 1 - \frac{n!×_{365}C_n}{365^n(365-n)!}\]

  • C は二項係数

\[_nC_k = \frac{n!}{k!(n-k)!}\]

  • R では choose(n, k)関数を使って計算できる　　
1 - ((factorial(50) * choose(365, 50)) / 365^50)
[1] 0.9703736
  • 50 人いるときに同じ誕生日のペアがいる確率は 0.97 (= 97%)
  • 人数ごとに同じ誕生日のペアがいる確率を求める式を関数化し、22 人と 23 人の時の確率を計算してみる
Birth_Expect <- function (n) {
  1 -  ((factorial(n) * choose(365, n)) / 365^n)
}
Birth_Expect(22)
[1] 0.4756953
Birth_Expect <- function (n) {
  1 - ((factorial(n) * choose(365, n)) / 365^n)
}
Birth_Expect(23)
[1] 0.5072972
  • 同じ誕生日のペアがいる確率が 50% を超える人数は 23 人
  • 「2.1.1 シミュレーション（試行回数 = 100）」で作成したデータフレーム (Prob_df) にこの理論値に Expect と名前を付けて追加する
Prob_df <- Prob_df |> 
  mutate(Expect = map_dbl(Students,
                          ~Birth_Expect(.x)))
DT::datatable(Prob_df)
  • pivot_longer()関数を使って tidy なデータに変換する
Prob_tidy <- Prob_df |> 
  pivot_longer(cols = Probs:Expect,
               names_to = "Type",
               values_to = "Prob")
DT::datatable(Prob_tidy)
  • シミュレーションから得られた結果 (Probs) と理論値 (Expect) を折れ線グラフにしてみる　　
Prob_tidy %>%
    mutate(Type = ifelse(Type == "Probs", "シミュレーション", "理論値")) %>%
    ggplot() +
    geom_line(aes(x = Students, 
                  y = Prob * 100, 
                  color = Type), 
              size = 1) +
    labs(x     = "学生数",
         y     = "同じ誕生日の人が2人以上いる割合 (%)\n(試行回数 = 100)",
         color = "") +
    theme_minimal(base_family = "HiraKakuProN-W3") +
    theme(legend.position = "bottom")

  • 上の図ではシミュレーションでの試行回数が 100回
  • 試行回数を 1000 回に増やして図を描いてみる
n_student <- 100    # 学生数
n_trials  <- 1000 # 試行回数

# 結果を格納する空ベクトルを用意する
Result_vec <- rep(NA, n_trials)

# 反復処理
for (i in 1:n_trials) {
    Birth_vec     <- sample(1:365, 
                            n_student, 
                            replace = TRUE)
    Result_vec[i] <- any(duplicated(Birth_vec))
}
Prob_df <- tibble(Students = 2:100,
                  Probs    = NA)

for (i in 1:nrow(Prob_df)) {
    
    Result_vec <- rep(NA, n_trials)
    
    for (j in 1:n_trials) {
        Birth_vec     <- sample(1:365, Prob_df$Students[i], replace = TRUE)
        Result_vec[j] <- any(duplicated(Birth_vec))
    }
    
    Prob_df$Probs[i] <- mean(Result_vec)
}
Prob_df %>%
    mutate(Expect = map_dbl(Students, ~Birth_Expect(.x))) %>%
    pivot_longer(cols      = Probs:Expect,
                 names_to  = "Type",
                 values_to = "Prob") |> 
    mutate(Type = ifelse(Type == "Probs", "シミュレーション", "理論値")) %>%
    ggplot() +
    geom_line(aes(x = Students, y = Prob * 100, color = Type), size = 1) +
    labs(x     = "学生数",
         y     = "同じ誕生日の人が2人以上いる割合 (%)\n(試行回数 = 1000)",
         color = "") +
    theme_minimal(base_family = "HiraKakuProN-W3") +
    theme(legend.position = "bottom")

  • モンテカルロ・シミュレーションから得られた結果と理論上の確率がほぼ近似していることが分かる

3. 確率と面積

3.1 円周率 \(π\) を使わずに円の面積を計算する

  • 円周率 \(π\) を使わずに、半径 1cm の円の面積を計算してみる
  • Step 1: 円の周りにぴったり当てはまる正方形を考える
  • Step 2: 正方形の中にランダムに点を打ちまくる
  • Step 3: 「正方形の中」と「円の中」に入った点の数の比を比べる
    →　「正方形の面積」と「円の面積」の比がわかる
    →　円の面積がわかる
    →　円周率 (3.14) がわかる

モンテカルロ・シミュレーションを使って、円周率を求めてみる

  • 正方形は辺の長さが 2 なので、正方形の面積 = 2 × 2 = 4
  • 正方形内に点をランダムに 100 個打ちまくる
set.seed(2022)
pi_df <- tibble(x = runif(100, -1, 1),
                y = runif(100, -1, 1))
pi_df
# A tibble: 100 × 2
         x      y
     <dbl>  <dbl>
 1  0.632   0.910
 2  0.295  -0.594
 3 -0.759   0.681
 4  0.0876  0.622
 5 -0.631   0.148
 6  0.272  -0.585
 7 -0.851   0.784
 8 -0.916  -0.831
 9 -0.259   0.243
10  0.515   0.968
# … with 90 more rows
pi_df |> 
  ggplot() +
    geom_rect(aes(xmin = -1, 
                  ymin = -1, 
                  xmax = 1, 
                  ymax = 1),
              fill = "white", color = "black") +
    geom_point(aes(x = x, y = y)) +
    coord_fixed(ratio = 1) +
    theme_minimal()

  • この正方形に円を追加する
  • 円を描くときには ggforce パッケージの geom_circle() 幾何オブジェクトを使う
  • マッピングは円の原点 (x = 0 と y = 0)、円の半径 (r)
  • 原点は (0, 0) で半径 \(r = 1\) の円を重ねる
pi_df %>%
    ggplot() +
    geom_rect(aes(xmin = -1, 
                  ymin = -1, 
                  xmax = 1, 
                  ymax = 1),
              fill = "white", color = "black") +
    geom_circle(aes(x0 = 0, 
                    y0 = 0, 
                    r = 1)) +
    geom_point(aes(x = x, y = y)) +
    coord_fixed(ratio = 1) +
    theme_minimal()

  • 円の中と外で点の色分けする
  • そのためには各点が円内に入っているかどうかを判定する変数 in_circle を作る
  • 点 (x, y) が、原点が (x′, y′)、かつ半径 r の円内に入っている場合、次の式が成立する

\[(x-x')^2 + (y - y')^2 < r^2\]

  • ここでは原点が (0, 0) で半径が 1 なので次の式で表せる

\[x^2 + y^2 < 1^2\]

  • 変数 in_circle を追加して、打ち込まれた点がこの条件を満たしているかどうかを判別する
pi_df <- pi_df %>%
    mutate(in_circle = if_else(x^2 + y^2 < 1^2, "円内", "円外"))

3.2 円の内と外に点が入る確率＝「面積」

  • 「円外」と「内側」に入った点の色を別々に表示させる
  • 散布図レイヤー geom_point() 内に color = in_circle と指定し、color を in_circle 変数でマッピングする
pi_df %>%
    ggplot() +
    geom_rect(aes(xmin = -1, 
                  ymin = -1, 
                  xmax = 1, 
                  ymax = 1),
              fill = "white", 
              color = "black") +
    geom_point(aes(x = x, 
                   y = y, 
                   color = in_circle), # color を in_circle 変数でマッピング
               size = 2) +
    geom_circle(aes(x0 = 0, 
                    y0 = 0, 
                    r = 1)) +
    labs(x = "X", y = "Y", color = "") +
    coord_fixed(ratio = 1) +
    theme_void(base_size = 12) +
   theme_minimal(base_family = "HiraKakuProN-W3")

  • group_by() 関数を使って、円内の点と円外の点の個数を数えてみる
pi_df |> 
  group_by(in_circle) |> 
  summarise(N = n())
# A tibble: 2 × 2
  in_circle     N
  <chr>     <int>
1 円内         82
2 円外         18
  • 100 個のうち、円内には 82 個 (82%)、円外には 18 個 (18%) の点が入った
  • 点の位置は無作為に決まる
  • 一片の長さが 2 の正方形の面積 = 2 × 2 = 4 なので
  • 100個の点は面積 4 の正方形の中に全部入るが、円の中に入った点は 82 個だから
    → 正方形の中に入っている円の面積は \(4 × \frac{82}{100} = 3.28\)
  • 円の半径が 1 なので、この値が円周率 \(π\)
4 * 0.82
[1] 3.28
  • 点の数を増やすとより正確な近似値が得られる
  • 2000個の点から得られた円周率を計算してみる
set.seed(2022)
pi_df2 <- tibble(x = runif(2000, -1, 1),
                 y = runif(2000, -1, 1))

pi_df2 <- pi_df2 %>%
    mutate(in_circle = if_else(x^2 + y^2 < 1^2, "円内", "円外"))

pi_df2 %>%
    ggplot() +
    geom_rect(aes(xmin = -1, 
                  ymin = -1, 
                  xmax = 1, 
                  ymax = 1),
              fill = "white", color = "black") +
    geom_point(aes(x = x, 
                   y = y, 
                   color = in_circle,
                   alpha = 0.5), 
               size = 2) +
    geom_circle(aes(x0 = 0, 
                    y0 = 0, 
                    r = 1)) +
    labs(x = "X", y = "Y", color = "") +
    coord_fixed(ratio = 1) +
    theme_void(base_size = 12) +
  theme_minimal(base_family = "HiraKakuProN-W3")

pi_df2 %>%
    group_by(in_circle) %>%
    summarise(N = n())
# A tibble: 2 × 2
  in_circle     N
  <chr>     <int>
1 円内       1562
2 円外        438
  • 2000 個のうち、円内には 1562 個、円外には 438 個
  • 点の位置は無作為に決まる
  • 正方形の面積 = 4 → 点が円内に入る確率 = 4 × = 円の面積
  • 円の半径が 1 なので、この値が円周率 \(π\)
4 * 1562 / 2000
[1] 3.124
  • だいぶ、3.14に近づいた

3.3 モンテカルロ・シミュレーション（円周率）

  • 打ち込む点の数 (x 軸) を 2 から 2000 まで増やした時に得られる円周率 (y 軸) をシミュレーションで求めてみる
dots <- 2:2000   # 点の数を 2 から 2000 まで変化させる
n <- length(dots)
results <- rep(NA, n)

for (i in 1:n) {
  results[i] <- tibble(x = runif(dots[i], -1, 1),
                       y = runif(dots[i], -1, 1)) |> 
    mutate(in_circle = if_else(x^2 + y^2 < 1^2, "円内", "円外")) |> 
    group_by(in_circle) |> 
    summarise(N = n()) |> 
    filter(in_circle == "円内") |>  # 円の内側に入った点だけを選ぶ
    summarise(mean(N/dots[i]*4))  # 推定された円周率を計算する
  } 

pi <- unlist(results)       # リストのデータをベクトルに変換
df_pi <- tibble(dots, pi) # データフレーム df_pi を作る

df_pi %>%
  ggplot() +
  geom_line(aes(x = dots, 
                y = pi), 
            size = 0.2,
            color = "blue",
            alpha = 0.8) +
  geom_hline(yintercept = 3.14, # 3.14 に赤い点線を引く 
             color = "red", 
             linetype = 2) +
  labs(x = "正方形に打ち込んだ点の数",
       y = "推定された円周率") +
  ggtitle("確率を使ったシミュレーションで得られた円周率") +
  theme_minimal(base_family = "HiraKakuProN-W3")

  • 点の数が増えるにつれて、安定して 3.14 が得られる傾向が見られる

4. ブートストラップ法

ブートストラップ法 (bootstrapping) とは

  • モデルが複雑すぎて統計量の標準誤差がうまく計算できない場合に使われ、推定量の不確実性の近似値を得る方法
  • Efron (1979) が提案した方法で、データ分析において広く使われている
  • ブートストラップにおいて重要なこと：
  • 元のデータセットのサンプルサイズを n とした場合
  • 復元抽出を用いてサンプルサイズ n のデータセットを再度構築　→　平均値を計算
  • この手順を 5000回繰り返せば　→　5000個の平均値が得られる
  • 5000個の平均値の平均値　→　元のデータセットの平均値に近似
  • 5000個の平均値の標準偏差 →　元のデータセットの平均値の標準偏差(= 標準誤差)に近似
  • 本当にこのようになるか試してみる
  • 長さ 100 のベクトルを 2 つ作成する
  • この 2 つのベクトルは平均値が 2、5、標準偏差 1 の正規分布に従うとする

\[data1 〜 Normal(\mu = 2, \sigma = 1)\\ data2 〜 Normal(\mu = 5, \sigma = 1)\]

  • 2 つの標本の平均値の差分は \(2-5= 3\)
set.seed(2020)
data1 <- rnorm(100, 2, 1)
data2 <- rnorm(100, 5, 1)
df <- data.frame(data1, data2)

df_long <- df %>% 
    tidyr::pivot_longer(data1:data2,
                 names_to = "data",  
                 values_to = "score") 
df_long
# A tibble: 200 × 2
   data   score
   <chr>  <dbl>
 1 data1  2.38 
 2 data2  3.27 
 3 data1  2.30 
 4 data2  4.01 
 5 data1  0.902
 6 data2  4.41 
 7 data1  0.870
 8 data2  5.38 
 9 data1 -0.797
10 data2  5.75 
# … with 190 more rows
df_long %>% 
  ggplot() +
  geom_boxplot(aes(x = data, 
                   y = score))

差分の不確実性

  • ここでは平均値の差の検定 (t検定)をやってみる
  • ここでは標準偏差が同じ (\(\sigma = 1\)) 2 つの分布から得られた標本
    →　等分散を仮定する検定を行う (var.equal = TRUE)
ttest_result <- t.test(df_long$score[df_long$data == "data1"],
                       df_long$score[df_long$data == "data2"],
                       var.equal = TRUE) 
ttest_result

    Two Sample t-test

data:  df_long$score[df_long$data == "data1"] and df_long$score[df_long$data == "data2"]
t = -17.423, df = 198, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -3.089598 -2.461319
sample estimates:
mean of x mean of y 
 2.108892  4.884350
  • 平均値の差分の不確実性 (= 標準誤差) は ttest_result$stderr で抽出可能
ttest_result$stderr
[1] 0.1592985
  • 平均値の差分の不確実性 (= 標準誤差) は約 0.159
  • 標準誤差の値さえ分かれば、検定統計量も、信頼区間も、値も計算できる
    →　重要なのは標準誤差
  • 試行回数は 1 万のブートストラップ法で標準誤差の近似値を計算する

①結果を格納する長さ 1 万の空ベクトル result_vec を作成
②i の初期値を 1 と設定
③data1 から 100個 (= data1 の大きさ) の値を無作為抽出→sample1 に格納
④data2 から 100個 (= data2 の大きさ)の値を無作為抽出 → sample2 に格納
⑤result_vec の i 番目の位置に sample1 の平均値と sample2 の平均値の差分を格納
⑥i を 1 増加させる
⑦③〜⑥の手順を 1 万回繰り返す

n_trials   <- 10000
result_vec <- rep(NA, n_trials)

for (i in 1:n_trials) {
    sample1 <- sample(data1, 100, replace = TRUE)
    sample2 <- sample(data2, 100, replace = TRUE)
    
    result_vec[i] <- mean(sample1) - mean(sample2)
}
  • result_vec には平均値の差分が 1万個格納されている
    →これらの値の平均値をブートストラップ推定量　デルタ \(\delta\) (bootstrap estimate) と呼ぶ
# 最初のデータの平均値の差分
delta <- mean(data1) - mean(data2)
delta
[1] -2.775459
# ブートストラップ推定量 (平均値の差分の平均値)
boot_delta <- mean(result_vec)
boot_delta
[1] -2.772783
# ブートストラップ推定量のバイアス
boot_b <- delta - boot_delta
boot_b
[1] -0.002676098
  • 実際の平均値の差分 (\(\delta_0\)) は -2.775459
  • ブートストラップ推定量 (\(\delta^*\)) は -2.772783
-2.775459 - (-2.772783)
[1] -0.002676
  • その差 -0.002676→ブートストラップ推定量のバイアス (\(b\))
  • 我々に興味があるのは平均値の差分ではない
  • 平均値の差分はブートストラップ法を用いなくても普通に計算できる
  • ここで知りたいのは平均値の差分の不確実性、つまり標準誤差
  • ブートストラップ推定量の標準誤差 = result_vec の標準偏差
# ブートストラップ推定量の標準誤差
boot_se <- sd(result_vec) 
boot_se
[1] 0.1579615
  • ブートストラップ推定量の標準偏差 (\(se\)) は約 0.1579615
  • t 検定の結果から得られた標準誤差 0.1592985 と非常に近い

95%信頼区間

百分位数信頼区間 (Percentile CI) は result_Vec の左側の領域が 2.5%、右側の領域が 2.5%となる区間
- quantile() 関数で計算できる

# 95%信頼区間
quantile(result_vec, c(0.025, 0.975))
     2.5%     97.5% 
-3.082493 -2.456123
  • バイアスを補正しない Wald 信頼区間を計算する
  • \(\delta_0 ± 1.96\)
    (1.96は標準正規分布の累積密度分布において下側領域が0.975となる点です。
delta + qnorm(c(0.025, 0.975)) * boot_se
[1] -3.085058 -2.465860
  • バイアス修正Wald信頼区間
(delta - boot_b) + qnorm(c(0.025, 0.975)) * boot_se
[1] -3.082381 -2.463184
  • 結果のまとめ
t検定 ブートストラップ
平均値の差分 -2.775458 -2.775459
標準誤差 0.1592985 0.1579615
95%信頼区間 [-3.09, -2.46]
95%信頼区間（百分位数） [-3.08, -2.46]
95%信頼区間 (Wald) [-3.09, -2.47]
95%信頼区間（バイアス修正 Wald） [-3.08, -2.46] 
2.108892-4.884350
[1] -2.775458

5. Monty Hall Problem

【ゲームのルール】

  • 3つのドア (1, 2, 3) に（高級車、ヤギ、ヤギ）がランダムに入っている
  • プレイヤーはドアを1つ選ぶ（例えば 1 を選んだとする）
  • モンティは残りのドアのうち 1 つを開ける（例えば 2 を開けたとする）
  • そのドアにはヤギが入っている（2 にはヤギが入っている）
  • プレイヤーはドアを選びなおすことができる （つまり、3 を選び直すことができる）
  • Question: プレイヤーは最初に選んだドアに止まるべきか、それともドアを選び直すべきか？

意見 A : ドアを選び直しても高級車を得られる確率は同じだから、どちらでも良い

意見 B: 高級車を得られる確率が上がるから、ドアを選び直すべき

  • この立場をとる説明の一例はこちら

【R を使ったシミュレーション】

  • 高級車はドア 1 にある

  • プレイヤーは三つのドア（1，2，３）のうち無作為に一つのドアを選ぶ
    （モンティがヤギの入っているドアを 1 つ開ける）

  • プレイヤーは次の二つの選択ができる：

    ① 最初に選んだドアにとどまる
    ② 残ったドアに選択を変える

    • プレイヤーは 200 回無作為にドアを選び、①と②それぞれの場合に高級車が当たる確率を計算する
# シミュレーションを行う為に 200個の入れものを準備。
set.seed(1838-3-11) # 乱数の設定。大隈重信の誕生日。
trials <- 200　 

# 「ドア 1」「ドア 2」「ドア 3」を無作為に 200 回選ぶ
prize <- sample(1:3, trials, replace = TRUE) 
prize
  [1] 3 1 3 1 2 2 2 1 1 3 3 2 1 1 2 1 2 3 3 3 3 3 1 2 3 1 2 3 1 1 2 1 1 3 2 1 2
 [38] 2 3 3 3 3 1 2 1 3 1 1 1 3 1 3 3 1 2 2 3 2 3 1 3 3 1 2 1 1 3 1 2 1 1 2 3 2
 [75] 3 1 3 2 2 2 1 3 2 2 3 3 2 3 2 2 2 2 1 3 3 2 3 1 2 2 3 2 3 2 3 2 3 1 1 2 1
[112] 2 2 3 1 1 1 3 3 1 3 2 1 2 1 2 3 1 2 2 2 2 1 3 1 2 1 2 2 3 3 3 2 1 1 2 3 3
[149] 2 3 2 1 2 3 3 2 3 3 1 1 3 2 3 3 2 3 2 3 3 3 3 2 1 1 3 3 3 2 2 1 3 2 3 3 1
[186] 3 2 3 3 1 2 2 1 2 1 1 1 3 3 2
  • 「ドア 1」に高級車があると指定
stay <- prize == 1    
head(stay)
[1] FALSE  TRUE FALSE  TRUE FALSE FALSE
  • stay では 1 が TRUE、2 と 3 が FALSE と表示
  • 「ドア 2」と「ドア 3」に高級車はないと指定
move <- prize != 1  
head(move)
[1]  TRUE FALSE  TRUE FALSE  TRUE  TRUE
  • move では 1 が FALSE、2 と 3 が TRUE と表示
  • Monty Hall の提案を受けず「ドア 1」に止まった時に高級車が当たる割合 prob.stay と「ドア 1」から動いた時に当たる割合 prob.move の累積を 200 回分計算する
  • シミュレーション結果をプロットする

Monty Hall Problem シミュレーションのためのスクリプト 

# シミュレーションを行う為に 200個の入れものを準備。
set.seed(1838-3-11) # 乱数の設定。大隈重信の誕生日。
trials <- 200　 
prize <- sample(1:3, trials, replace = TRUE) # 「ドア 1」「ドア 2」「ドア 3」を無作為に 200 回選ぶ
prize
stay <- prize == 1 #「ドア 1」に高級車があると指定
head(stay)
move <- prize != 1 #「ドア 2」と「ドア 3」に高級車はないと指定
head(move)
prob.stay <- prob.move <- rep(NA, trials) 
for (i in 1:trials) {
prob.stay[i] <- sum(stay[1:i]) / i     # 「ドア 1」に止まった時に高級車が当たる割合
prob.move[i] <- sum(move[1:i]) / i 　　# 「ドア 1」から動いた時に高級車が当たる割合

plot(1:trials, prob.move, ylim = c(0, 1),  #シミュレーション結果をプロット
     type = "l", col = "red", 
     xlab = "number of trials",
     ylab = "prob. of getting a car",
     main = "Simulation of Monty Hall Problem")
 lines(1:trials, prob.stay, type = "l", col = "blue")
 abline(h = c(0.33, 0.66), lty = "dotted") # 勝つ確率が 33% (1/3)と66% (2/3)に点線を引く
 legend("topright", lty = 1, col = c("red", "blue"),
       legend = c("change", "stay"))}

・実行方法
・R 上で「ファイル」=> 「新規文書」
=> 上のコード（Monty Hall Problem シミュレーションのためのスクリプト）をペースト
=> Return (Enter) を押す

注意：RMarkdown上で実行できないこともないが、動画アニメーションではなく「紙芝居」のような動きをするので、R 上で実行するのがよい

参考文献