問4はベイズ推定を含んだ比較的シンプルな問題でした。統計数理で出てもよさそうな内容なのでざっと解いておくと良さそうです。
(1)
まずは普通に二項分布の最尤推定を行う問題です。Y=yの下での尤度関数は以下になります。
\(L(\pi|Y=y)=_nC_y\pi^y(1-\pi)^{n-y}\)
ここから対数尤度関数をとって、微分しますと
\(\frac{\partial}{\partial\pi}logL=\frac{y}{\pi}-\frac{n-y}{1-\pi}\)
これ＝０とおいて最尤推定値を求めると
\(\hat\pi=\frac{y}{n}\)
となります。
(2)
さてここからはベイズです。事後分布を求めていく問題。求めたい事後分布は以下の式で表されます。
\(f(\pi|y)=\frac{P(y|\pi)f(\pi)}{P(y)}\)
ここで分子の\(P(y|\pi)\)は先ほどの尤度関数と同じ式で表すことができ、分母については以下のようにπで周辺化することでP(y)を求めます。
\(P(y)\\=\int_0^1{_nC_y}\pi^y(1-\pi)^{n-y}d\pi\)
さてこれを分子と比較すると共通しているコンビネーションの部分は消えますので
結局は事後分布が以下の式になります。
\(f(\pi|y)\\=\frac{\pi^y(1-\pi)^{n-y}}{\int_0^1\pi^y(1-\pi)^{n-y}d\pi}\\=\frac{\pi^y(1-\pi)^{n-y}}{B(y+1, n-y+1)}\)
分母が実はβ関数になっていますので、上記のようなβ分布に変形できます。
実はこれは当然といえば当然の結果です。そもそも一様分布はベータ分布の特殊形で
beta(1,1)の場合でした。
そして、ベータ分布は共役事前分布（事後分布になっても分布の形は同じとなる）であるので
再びベータ分布となります。
(3)
最頻値は事後分布を微分し、０になるところが頂点となりますのでそれを確認します。
\(f'(\pi|y)\\=y\pi^{y-1}(1-\pi)^{n-\pi}-(n-y)\pi^y(1-\pi)^{n-y-1}\\=0\)
とおいて計算すると
\(\pi=\frac{y}{n}\)
となり、最尤推定値と合致することが分かりました。
(4)
事後分布の期待値を求める問題です。ベータ分布の期待値は一般にB(a,b)に対してa/a+bとなりますので、答えはB(y+1,n-y+1)に対して
\(E[\pi|y]=\frac{y+1}{n+2}\)
となります。なお、これは最尤推定量y/nと事前分布の期待値1/2の重み付き平均をとったものとなっています。（最尤推定量の重みがサンプルサイズなのはわかるのですが、事前分布の期待値の２はどこから来たのか？が疑問です）
ついでと言っては何ですがベータ分布の期待値の導出も一般に従って書いておきます。
＜補足＞ベータ分布の期待値
前提としてベータ関数とガンマ関数の以下の関係を使います。
\(B(a,b)=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a,b)}\)
なお、導出は以下の記事で紹介しています。久々に見ましたが結構導出は大変。
さて、ここでbeta(a,b)に対して期待値を計算すると(a>0, b>0, 0<x<1)
\(E[x]\\=\int_0^1\frac{x^{a}(1-x)^{b-1}}{B(a,b)}dx\\=\frac{B(a+1,b)}{B(a,b)}\int_0^1\frac{x^{a}(1-x)^{b-1}}{B(a+1,b)}dx\\=\frac{B(a+1,b)}{B(a,b)}\\=\frac{\Gamma(a+1)\Gamma(b)}{\Gamma(a+b+1)}\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\\=\frac{a}{a+b}\)
となります。
(5)
ｙの周辺分布を求める問題です。πの全範囲で積分してやりましょう。
\(\int_0^1P(y|\pi)d\pi\\=_nC_y\int_0^1\pi^y(1-\pi)^{n-y}d\pi\\=_nC_yB(y+1, n-y+1)\)
ここでｙは自然数ですのでガンマ関数に変換後、階乗に変換できることが分かります。
\(_nC_yB(y+1, n-y+1)\\=\frac{n!}{(n-y!)y!}\frac{y!(n-y!)}{(n+1)!}\\=\frac{1}{n+1}\)
y=0,1,…,nをとる離散一様分布であることが分かりました。
(6)
実際に事後分布を使って確率計算をしてみようという問題です。当然ｙが大きい方がπが1/2以上となる事後確率は上がりますので、数値の大きめな方を狙って計算してみるのが何となく良いでしょう。
計算方法ですが、πが1/2以上である確率は事後分布より以下の式で表されます。
\(P(\pi\geq\frac{1}{2}|y)\\=\frac{1}{B(y+1,6-y)}\int_{\frac{1}{2}}^1\pi^y(1-\pi)^{5-y}\)
あとはひたすら数値を当てはめてみますと
y=4のとき0.89, y=3のとき0.65となりましたので
最小値は４となります。
統計学をどんどん自分で学んで深めたい、という方へのおすすめ書籍をこちらのページの下部にまとめております。初心者向けのものから応用まで幅広く読んでいますので参考にどうぞ。
コメントを残す