実習7.3-4 複数の山はムリ、h(θ) は分母の近似

「実習7.1-2 ランダムウォークを観察」からの続き。


実習7.3 コインの状態が「不正なし、表か裏に偏っている」という事前確率を、 [cos(4πθ) + 1]² （θ：表の出る確率）とする。観測データは 12 回の試行で 8 回の表。

このコインの事前確率は、実習6.2 と同じだが、今回の実装方法は cos 関数を利用。ベイズルールより

これが事後確率確率。

ところが、この事後確率と事前確率 [cos(4πθ) + 1]² が同じ種類の分布ではない。それに、事前確率 [cos(4πθ) + 1]² と尤度 θ⁸ × (1 - θ)⁴ が conjugate（結合）できるとは思えない。

以下、本書提供のスクリプト BernGrid.R と BernMetropolisTemplate.R の実行例を示す。

BernGrid.R

> binwidth<-1/1000
> thetagrid<-seq(from=binwidth/2,to<-1-binwidth/2,by=binwidth)
> relprob<-(cos(4*pi*thetagrid)+1)²
> prior<- relprob / sum(relprob)
> datavec<-c(rep(1,8),rep(0,4))
> BernGrid(Theta = thetagrid,pTheta = prior,Data = datavec)

BernMetropolisTemplate.R

スクリプトの修正点は以下の通り

myData = c( 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0 )
trajLength = 111111
prior = function( theta ) {
  prior <- (cos(4*pi*theta)+1)^2
  prior[theta>1 | theta<0] <- 0
  return( prior )
}


どちらの実行結果も「2 つの山」を持つ分布だが、HDI の前提は

unimodal distribution （山が1つの分布）

よって、この事前確率（ [cos(4πθ) + 1]²）を分析する「単純」な方法はない。

「ではどうするのか？」は、今のところ私も分かりません。今後の課題です。

実習7.4 BernMetropolisTemplate.R の h(θ) 変更後の weighted evidence を見る。

「Metropolis Algorithm の実装例」で解説したように、ステップ i 毎の weighted evidence は左のようになる。今回の h(θ) は ベータ関数で、サンプリングした事後確率の平均と標準偏差が、h(θ) の平均と標準偏差となる。

BernMetropolisTemplate.R の末尾に以下を追加して実行

　X11(); plot(wtdEvid,type="l")

これは、9,999 ステップ（X 軸の Index）毎の weighted evidence をプロットしたもの。wtdEvid の内訳は

> summary(wtdEvid)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
   4003    5424    5521    5463    5561    5573

また、事後確率に基づいて算出した a, b は
> a; b
[1] 12.38976
[1] 4.154175

次に、この a, b に事後確率とは無関係の値に変更して実行する。

左は a = 1, b = 1 、右は a = 10, b = 10 。両方とも Y軸の値が、変更前とあまりにも大きい。

h(θ) が、事後確率を近似する必要性がわかる。


「実習7.5 モデル 3 つの比較」に続く。