BDA3の3.10 Exercisesの1を解く。 本文はBayesian Data Analysis。 解答参考はSolutions to some exercises from Bayesian Data Analysis。
観測値(y1,…,yJ)は多項分布Multinomial(n,(θ1,…,θJ))に従う。 多項分布Multinomial(10,(0.6,0.2,0.2))は赤6個、青2個、黄2個のボールが10個入った箱から ランダムに取り出した時に、それぞれが取り出される回数の分布である。
using Distributions, StatsPlots, Random
Random.seed!(1)
rand(Multinomial(10, [0.6,0.2,0.2]))
3-element Vector{Int64}:
7
0
3
シード値1で赤7回、青0回、黄3回が取り出された。
今、θの事前分布をディリクレ分布Dirichlet(a1,…,aJ)とする。 ディリクレ分布の確率密度関数は
p(θ∣a)∝j=1∏Jθaj−1 である。多項分布の確率密度関数は
p(y∣θ)∝j=1∏Jθjyj であり、ディリクレ分布は多項分布の共役事前分布となっている。 θの事後分布はDirichlet(y1+a1,…,yJ+aJ)となる。
θをθ1、θ2、θr=1−θ1−θ2で分ける。θrはθ1、θ2以外の全てを足したものである。 確率密度関数は
p(θ1,θ2∣y)∝θ1y1+a1−1θ2y2+a2−1(1−θ1−θ2)yr+ar−1 となる。
α=θ1+θ2θ1 β=θ1+θ2 とおく。
α=θ1+θ2θ1(θ1+θ2)α=θ1θ1=αβ β=θ1+θ2β=αβ+θ2θ2=(1−α)β と変換できる。
(α,β)=(θ1+θ2θ1,θ1+θ2) 変換のヤコビアンJを計算する。
J=(∂θ1∂α∂θ1∂β∂θ2∂α∂θ2∂β) であるので、それぞれ計算し、
J=((θ1+θ2)2θ21(θ1+θ2)2−θ11) ∣J∣=(θ1+θ2)2θ2−(θ1+θ2)2−θ1=θ1+θ21=β1 となる。
α、βの周辺事後分布を求めるために、(3)をα、βに変換する。
p(α,β∣y)∝∣J∣1(αβ)y1+a1−1((1−α)β)y2+a2−1(1−β)yr+ar−1=βαy1+a1−1(1−α)y2+a2−1βy1+a1−1+y2+a2−1(1−β)yr+ar−1=αy1+a1−1(1−α)y2+a2−1βy1+a1+y2+a2−1(1−β)yr+ar−1 α部分とβ部分に分けることができ、αはBeta(y1+a1,y2+a2)であり、 βはBeta(y1+a1+y2+a2,yr+ar)である。
最初は多項分布だったが、θ1、θ2の2つの確率に注目するためにαを導入している。 その結果、周辺事後分布はベータ分布で表すことができ、y3,…,yJにも影響を受けないことが分かった。