The One with ...

思いついたことや作業メモなどを公開しています

Stanによる推定例:ベルヌーイ分布のパラメータ

連勝記録から勝率を推定する

ある棋士が29連勝した,というデータから,勝利確率の事後分布をMCMCで計算してみる.統計モデルの仮定は次の通りである.

  • 毎回の勝ち負けは独立にベルヌーイ分布に従っている

  • ベルヌーイ分布のパラメータpの事前分布は範囲[0,1]の一様分布とする

StanとRのコードは以下の通り(コードを書くうえで松浦健太郎さんの『StanとRでベイズ統計モデリング』を参考にさせていただきました)

パラメータpの事後分布のプロット.

f:id:hamada7418:20170628100158p:plain

ベルヌーイ分布にしたがう確率変数の例を日常生活の中から見つけてみよう.

ちなみに勝率0.9で29連勝する確率は

0.929=0.04710129

である.

解析的な事後分布の導出

上述の計算機による推定と同様の結果を,手計算によって導出してみよう.

試行数1のベルヌーイ分布の確率関数は

{ \displaystyle
p ( Y | \theta) = \theta^{Y} (1-\theta)^{1-Y}
}

だから,n回試行の同時確率は

{ \displaystyle
\prod_{i=1}^{n} \theta^{Y_{i}} (1-\theta)^{1-Y_{i}}
}

である.この同時関数をパラメータθの関数として見たものが尤度関数{p(y|\theta) }である. θの事前分布としてベータ分布

{ \displaystyle
B(a,b)
}

を仮定する.さて事後分布は

{ \displaystyle
p(\theta|y) \propto p(\theta)p(y|\theta)
}

{ \displaystyle \quad \propto \theta^{a-1} (1-\theta)^{b-1} \prod_{i=1}^{n} \theta^{Y_{i}} (1-\theta)^{1-Y_{i}}
 }

であった.ここで積の前半部分はベータ分布のカーネル 

{ \displaystyle \theta^{a-1} (1-\theta)^{b-1} } 

であり,後半部分はベルヌーイ分布の尤度

{ \displaystyle \prod_{i=1}^{n} \theta^{Y_{i}} (1-\theta)^{1-Y_{i}} }

である.その計算結果は

{ \displaystyle
\theta^{\alpha-1} (1-\theta)^{\beta-1}  }

ただし { \displaystyle \alpha = a + \sum_{i=1}^{n}y_{i}, \quad \beta = b + n - \sum_{i=1}^{n}y_{i} } 

となる.これはベータ分布のカーネルである.つまり事前分布としてベータ分布を指定すると,ベルヌーイ分布のパラメータの事後分布はベータ分布

{ \displaystyle
\theta | y \sim B( a + \sum_{i=1}^{n}y_{i},  b + n - \sum_{i=1}^{n}y_{i})
 }

となる.よって事前分布としてa=1, b=1と仮定すればθの事後分布は

{ \displaystyle
B(1+29, 1+29-29) = B(30, 1)
 }

である.このベータ分布をθの事後確率分布としてplotとしてみよう.まず事前分布がB(1,1)の場合.

f:id:hamada7418:20170713000704p:plain

次に事前分布がB(2, 2)の場合.

f:id:hamada7418:20170713000707p:plain

Stanで推定した結果に対応した確率密度関数のplotを得た*1

*1:厳密に言うと,Stanの結果のplotはMCMCサンプルのヒストグラムなので,確率密度関数のplotではない.