Juliaの世界

BDA3の3.10 Exercisesの1を解く。 本文はBayesian Data Analysis。 解答参考はSolutions to some exercises from Bayesian Data Analysis

  1. Binomial and multinomial models

観測値(y1,,yJ)(y_1,\dots,y_J)は多項分布Multinomial(n,(θ1,,θJ))\mathrm{Multinomial}(n,(θ_1,\dots,θ_J))に従う。 多項分布Multinomial(10,(0.6,0.2,0.2))\mathrm{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)\mathrm{Dirichlet}(a_1,\dots,a_J)とする。 ディリクレ分布の確率密度関数は

p(θa)j=1Jθaj1 p(θ|a)∝\prod^J_{j=1}θ^{a_j−1}

である。多項分布の確率密度関数は

p(yθ)j=1Jθjyj p(y|θ) ∝ \prod^J_{j=1}θ^{y_j}_j

であり、ディリクレ分布は多項分布の共役事前分布となっている。 θの事後分布はDirichlet(y1+a1,,yJ+aJ)\mathrm{Dirichlet}(y_1+a_1,\dots,y_J+a_J)となる。

θの事後分布

θをθ1θ_1θ2θ_2θr=1θ1θ2θ_r=1-θ_1-θ_2で分ける。θrθ_rθ1θ_1θ2θ_2以外の全てを足したものである。 確率密度関数は

p(θ1,θ2y)θ1y1+a11θ2y2+a21(1θ1θ2)yr+ar1 p(θ_1,θ_2|y)∝θ_1^{y_1+a_1-1}θ_2^{y_2+a_2-1}(1-θ_1-θ_2)^{y_r+a_r-1}

となる。

(a) αの周辺事後分布

α=θ1θ1+θ2 α=\frac{θ_1}{θ_1+θ_2} β=θ1+θ2 β=θ_1+θ_2

とおく。

α=θ1θ1+θ2(θ1+θ2)α=θ1θ1=αβ α=\frac{θ_1}{θ_1+θ_2} \\ (θ_1+θ_2)α=θ_1 \\ θ_1=αβ β=θ1+θ2β=αβ+θ2θ2=(1α)β β=θ_1+θ_2 \\ β=αβ+θ_2 \\ θ_2=(1-α)β

と変換できる。

(α,β)=(θ1θ1+θ2,θ1+θ2) (α, β)=(\frac{θ_1}{θ_1+θ_2},θ_1+θ_2)

変換のヤコビアンJを計算する。

J=(αθ1αθ2βθ1βθ2) J=\begin{pmatrix} \frac{∂α}{∂θ_1} & \frac{∂α}{∂θ_2} \\ \frac{∂β}{∂θ_1} & \frac{∂β}{∂θ_2} \\ \end{pmatrix}

であるので、それぞれ計算し、

J=(θ2(θ1+θ2)2θ1(θ1+θ2)211) J=\begin{pmatrix} \frac{θ_2}{(θ_1+θ_2)^2} & \frac{-θ_1}{(θ_1+θ_2)^2} \\ 1 & 1 \\ \end{pmatrix} J=θ2(θ1+θ2)2θ1(θ1+θ2)2=1θ1+θ2=1β |J|=\frac{θ_2}{(θ_1+θ_2)^2}-\frac{-θ_1}{(θ_1+θ_2)^2}=\frac{1}{θ_1+θ_2}=\frac{1}{β}

となる。

α、βの周辺事後分布を求めるために、(3)をα、βに変換する。

p(α,βy)1J(αβ)y1+a11((1α)β)y2+a21(1β)yr+ar1=βαy1+a11(1α)y2+a21βy1+a11+y2+a21(1β)yr+ar1=αy1+a11(1α)y2+a21βy1+a1+y2+a21(1β)yr+ar1 p(α,β|y)∝\frac{1}{|J|}(αβ)^{y_1+a_1-1}((1-α)β)^{y_2+a_2-1}(1-β)^{y_r+a_r-1} \\ =βα^{y_1+a_1-1}(1-α)^{y_2+a_2-1}β^{y_1+a_1-1+y_2+a_2-1}(1-β)^{y_r+a_r-1} \\ =α^{y_1+a_1-1}(1-α)^{y_2+a_2-1}β^{y_1+a_1+y_2+a_2-1}(1-β)^{y_r+a_r-1}

α部分とβ部分に分けることができ、αはBeta(y1+a1,y2+a2)\mathrm{Beta}(y_1+a_1,y_2+a_2)であり、 βはBeta(y1+a1+y2+a2,yr+ar)\mathrm{Beta}(y_1+a_1+y_2+a_2,y_r+a_r)である。

(b) 結果の解釈

最初は多項分布だったが、θ1θ_1θ2θ_2の2つの確率に注目するためにαを導入している。 その結果、周辺事後分布はベータ分布で表すことができ、y3,,yJy_3,\dots,y_Jにも影響を受けないことが分かった。