【Rで検証編】「未完のゲーム」を終わらせよう〜「数学ガールの秘密ノート／確率の冒険」研究問題5-2
カバー写真はUnsplashのKarolina Kołodziejczakが撮影した写真。ダジャレですけど。
数学ガール
結城浩先生の「数学ガールの秘密ノート 確率の冒険」を読んで、巻末の「研究問題」についてあれこれ考えています。
前回は、偏ったコインを用いた「未完のゲーム」について考えました。これが果たして正しい結果を出してくれているのかについて、Rを用いて検証していきます。
研究問題5-2
研究問題5-2は、こういう問題です。
第5章〈未完のゲーム〉では、フェアなコインを投げてゲームを行いました。もしも、偏ったコイン（表が出る確率が1/2ではないコイン）を用いた場合は、どのような答えになるでしょうか。
そして、次のような答えを導き出しました。
$$
\displaystyle P(a,b)=\sum_{k=0}^{b−1}{n \choose k} p^{n−k}(1−p)^k,\quad (n=a+b−1)
$$
公正なコインを用いた場合は、式中のpが0.5の場合にあたり、$${p^{n-k}(1-p)^k}$$が常に$${\frac{1}{2^n}}$$になるため、次の式のように簡略化できます。
$$
\displaystyle \frac{1}{2^n} \sum_{k=0}^{b-1} {n \choose k}
$$
では、この式をRのコマンドに翻訳していきましょう。
Rによる検証
式をそのまま翻訳してみる
上の式（p=0.5の場合）を、そのまま関数に翻訳してみました。「ufg」は、unfinished gameの意味です。
# 解決した数式通りに書いてみる
ufg <- function(a,b) {
n <- a+b-1
p <- sum(choose(n, 0:(b-1))) / (2^n)
return (p)
}
# 実行例
ufg(2,2) # [1] 0.5
ufg(3,2) # [1] 0.3125まず、解決した数式通りに関数を組んでみます。関数choose(n, k)が二項係数を計算する関数です。
上のスクリプトでは、引数のkの部分に〈0 : (b-1)〉というベクトルが入れてあります。b=2のとき、〈0 : (b-1)〉は [0, 1] というベクトルと解釈されますから、choose(n, 0:(b-1))は、choose(n, 0)とchoose(n, 1)の結果をベクトルで返してくるのです。
choose関数の結果をsum関数で合計し、それを$${2^n}$$で割った値が求める確率です。
実行すると、実行例のような結果が返ってきます。ufg(2, 2)=0.5 は、本文でも確かめられていましたが、ufg(3, 2)=0.3125 についても、確かめておきましょう。
P(3,2)=0.3125の検証（p=0.5のとき）
Aが勝つのは、k=0,1の場合です。
k=0のとき：$${\frac{1}{2^4}{4 \choose 0} = \frac{1}{2^4}=0.0625}$$
k=1のとき：$${\frac{1}{2^4}{4 \choose 1} = \frac{4}{2^4}=0.25}$$
合計すると、$${0.0625+0.25=0.3125}$$となります。上の関数は正しい答えを返してくれていることがわかります。
二項分布関数を使って書いてみる
ところで、せっかくRを使っているので、二項分布関数を使ってより簡単に書き直してみましょう。
# 二項分布関数を使って書いてみる
ufg1 <- function(a, b, prob) {
return (sum(dbinom(0:(b-1), a+b-1, 1-prob)))
}
# 実行例
ufg1(1,1,0.5) # [1] 0.5
ufg1(2,2,0.5) # [1] 0.5
ufg1(3,2,0.5) # [1] 0.3125二項分布の確率を返す関数 dbinom(x, size, prob) を使います。この関数の引数は、probは表の出る確率で0.5、sizeがコイン投げの回数で式中ではn、xはそのうち表の出る回数で、式中ではn-kと解釈できます。したがって、この関数をそのまま使うと、Bが勝つ確率を返してしまいます。そこで、下のスクリプトでは、probを1-probと書き換えることで、Aが勝つ確率を返すようにしています。
実行例にあるように、正しい確率を返していることがわかります。これで、確率を変えて実行すれば、研究問題5-2に答えることができます。
確率を変えてみる
では、実際に実行してみましょう。
ufg1(1,1,0.6) # [1] 0.6
ufg1(2,2,0.6) # [1] 0.648
ufg1(3,2,0.6) # [1] 0.4752まず、P(1, 1)は、表が出る確率と同じになるはずですから、P(1, 1)=0.6 となります。では、P(2, 2) を確かめてみましょう。
P(2,2)の検証（p=0.6のとき）
Aが勝つのは、k=0,1の場合です。
k=0のとき：$${{3\choose0}\times0.6^3=0.216}$$
k=1のとき：$${{3\choose1}\times0.6^2\times 0.4=0.432}$$
したがって、これを合計して$${0.216+0.432=0.648}$$ となります。正しい答えが返っています。
P(3,2)の検証（p=0.6のとき）
Aが勝つのは、k=0,1の場合です。
k=0のとき：$${{4\choose0}\times0.6^4=0.1296}$$
k=1のとき：$${{4\choose1}\times0.6^3\times 0.4=0.3456}$$
したがって、これを合計して$${0.1296+0.3456=0.4752}$$ となります。正しい答えが返っています。
次回はアニメーション版
さて、次回はアニメーション版について書きます。でも、あまり面白くありません。なぜかというと、
この設定ならこの確率になるはず、という「数式による正解」が直感的にわからない。（計算しなくてはいけない）
計算した答えが、たとえば0.4752のように中途半端な値であるため、シミュレーションの結果がその値に収束していっているのか、とてもわかりにくい。（0.5より少し小さいことくらいしか結局はわからない）
というのが主な理由です。先に言い訳をしておく戦略でした。