2015年10月10日土曜日

ベイズの基礎:課題16-3 コインが表の確率分布

課題16-2 何割打者か? からの続き。


Activity 16-3: Spinning a penny

課題:コインを「立てて」回転させて裏か表の結果を記録する。この試行を5回行い「お表が出る確率」を求める。

予測モデルは 0, 0.1, 0.2, ..., 0.9, 1.0 の 11 個。の場合は「表は出ない」、1.0 の場合は「表しか出ない」という意味。そうなると完全に「不正コイン」だが、0.1 や 0.9 でも十分「不正コイン」。

事前確率は、そのコインに対する予備知識がないので、11 個のモデルは平等に有り得るとして全て 1/11 = 0.091 とする。

尤度の算出方法は、例えば {表, 裏, 表, 裏, 表} の出る確率は「モデル 0.4」の場合

尤度 for モデル0.4 = 0.4 × (1 - 0.4) × 0.4 × (1 - 0.4) × 0.4

次は「表の出方」は全6通り、その全てで事後確率を求めたのが以下の表。
「表が2回」の事後確率を「確率分布」でグラフ化。

左が「事前確率」1/11 = 0.091 の一様分布。これが観察データにより分布が変わったのが右。表が2回なので、若干左(低い予測モデル)の方に偏っていることが分かる。コインに不正がなければ、試行回数を増やすごとに「0.5モデル」をピークにした正規分布に近づく。

y = c(0.0000,0.0437,0.1229,0.1852,0.2074,0.1875,0.1383,0.0794,0.0307,0.0049,0.0000)
par(mfrow=c(1,2))
plot(seq(0,10), rep(0.091,11), type="h", main="PRIOR", xlab="", ylab="", ylim=c(0, 0.22))
plot(seq(0,10), y, type = "h", main="2 Heads in 5 Spins", xlab = "", ylab="", ylim=c(0, 0.22))

左は、確率の高いモデルから順に並べて、確率の累積値を求めたもの。

この表から「50% estimate(50%見積り)区間」や「90% estimate(90%見積り)区間」が分かる。

58%見積り区間:{0.4, 0.5, 0.3} つまり
92%見積り区間:{0.4, 0.5, 0.3, 0.6, 0.2, 0.7}

これらの区間で、最大と最小のモデルから

 [0.3, 0.5] は58% probability interval(58%の確率区間)
 [0.2, 0.7] は92% probability interval(92%の確率区間)

少し分かりにくいかもしれない。上記の事後確率のグラフと、確率の累積値を求めた表から理解できると思います。

左は累積値をグラフ化したもの。50%を青い点線、90%を赤い点線で記した。

plot(cumsum(sort(y, decreasing = T)), type = "o")
abline(h=0.9, lty=2, col="red")
abline(h=0.5, lty=2, col="blue")


課題16-4 夫が妻より年上な確率 に続く。

0 件のコメント:

コメントを投稿