大数の法則を確かめよう〜別解？というか、むしろこっちが正解？

2024年9月27日 14:01

前回投稿した次の記事で、「もしかしたら題意を読み違えていたかも？」と思ったので、別解を書いておきます。別解と言っていますが、こっちのほうがより題意に沿ったものかもしれません。

問題の確認

問題はこうでした。

実際にコインを投げて、表が出るか裏が出るかを調べましょう。コインをM回まで投げた時点で表が出た回数(m)を集計し、横軸がMで縦軸が相対度数m/Mとなるグラフを描いてみてください。

「数学ガール確率の冒険」より（一部要約、改変）

「コインをM回まで投げた時点で表が出た回数（ｍ）を集計し」とあるので、「10回まで投げたぞ。ええと、表は4回で、相対度数は0.4だ」という記録があり、そのあと続けてまた10回投げ、「合計で20回投げたぞ。ええと、表が11回で、相対度数は0.55だ」という記録がある、というのが、題意なのではないか。

何が違うか

前回投稿したものは、「まず10回投げて、表が出る割合を計算する」。次に、その記録はそのまま置いておいて、「新しく20回投げて、表が出る割合を計算する。」さらに、その記録とは別に、「新しく30回投げて、・・・」という方法でグラフを書いていました。
しかし、上に書いたのが正しい題意だとして、それにそってプログラムを書くなら、「まず10回まで投げて、相対度数を計算する」、「続けてもう10回投げて、合計で20回までの相対度数を計算する」、というようにするのが良いでしょう。

書き直したコード

というわけで、コードを書き直しました。ついでなので（何のついでかよくわからないが）、5人いっしょに実験したら楽しかろうと思いついて（楽しいのか、本当に？）、5本の折れ線グラフを重ね描きするコードに変えてみました。

# 別解
# フェアなコイン：表が出る確率を0.5とする
prob <- 0.5
set.seed(1234)
# サイコロ投げ1回から5000回まで試す
trial_max <- 5000
# 友達といっしょに実験します
person <- 5
# データを記録するデータフレーム
dat <- data.frame(id=1:trial_max)
for (i in 1:person) {
  cn1 <- paste("Dice", i, sep="_")
  cn2 <- paste("Ratio", i, sep="_")
  dat[, cn1] <- NA
  dat[, cn2] <- NA
}
# サイコロ投げ実験開始
for (t in 1:trial_max) {
  for (i in 1:person) {
    dat[t, i*2] <- rbinom(1, 1, prob)
    dat[t, i*2+1] <- mean(dat[1:t, i*2])
  }
}
# データをまとめる
cn <- seq(1, (person+1)*2, 2)
db <- subset(dat[, cn])
# グラフ
matplot(db$id, db[, -1], type="l", ylim=c(0,1), 
        xlab="Trials", 
        ylab="Proportion", 
        main=paste("Probability =",prob))
abline(h=prob, col="red", lty=2)
# legend("topright", legend=1:person, lty=1:person, col=1:person)

面倒なのでコードの説明はしません。（おい）
途中、データを記録するためのデータフレームを作っています。最後に matplot 関数で折れ線グラフを描くときに、このほうが便利なので、わざわざデータフレームにしていますが、もっとスマートな方法がきっとあるんだろうと思います。
では、結果です。

最初のうち、1000回くらいまでは、表の出る割合が0.45から0.55くらいまでの幅をもっていますが、それ以降範囲がせばまって、2000回を超える頃から、ほとんど変動しなくなっていきます。

今回のコードでは、乱数シードを固定していますから（set.seed(1234)の行がそれです）、同じシードを使えば同じ結果が得られます。別のシードもいくつか試してみましたが、おおむね似たような結果になりますが、0.5への近づき方は、結構ばらばらでした。

今回は以上です。Javascriptを使ったシミュレーションは、今回のコードをもとにして行う予定です。

カバー画像は前回と同じで

カバー写真はUnsplashのAndre Taissinが撮影した写真でした。