Juliaの世界

BDA3の2.11 Exercisesの17を解く。 リンクはBayesian Data Analysis

  1. Posterior intervals

最高事後信用区間(highest posterior interval)は中央事後信用区間(central posterior interval)と違い、変換に関し不変でない。 nνσ2χn2\frac{nν}{σ^2} \sim χ^2_nとし、σは無情報事前分布p(σ)σ1p(σ)∝σ^{-1}を持つとする。

using Distributions, StatsPlots
plot(Chisq(10),yrotation=90,label="Chisq(10)")
savefig(joinpath(@OUTPUT, "Chisq10.svg"))
nothing

χn2χ^2_nは左右対称でない。

(a) σ2σ^2の事前分布

u=σ2u=σ^2とおくとσ=uσ=\sqrt{u}である。

p(σ2)=p(u) p(σ^2)=p(u)

確率密度関数は1対1の変数の変換φ=h(θ)φ=h(θ)に関し

p(φ)=p(θ)dθdφ=p(θ)h(θ)1 p(φ)=p(θ)|\frac{dθ}{dφ}|=p(θ)|h'(θ)|^{-1}

が成り立つ。今、θ=uθ=\sqrt{u}φ=h(θ)=θ2=uφ=h(θ)=θ^2=uとすると、

p(u)=p(u)dudu=p(u)12u12 p(u)=p(\sqrt{u})|\frac{d\sqrt{u}}{du}|=p(\sqrt{u})\frac{1}{2}u^{-\frac{1}{2}}

u=σ2u=σ^2σ=uσ=\sqrt{u}なので

p(σ2)=12σ1p(σ) p(σ^2)=\frac{1}{2}σ^{-1}p(σ)

σは事前分布として、p(σ)σ1p(σ)∝σ^{-1}であると仮定しているので、

p(σ2)σ1σ1=σ2 p(σ^2)∝σ^{-1}σ^{-1}=σ^{-2}

(b) 最高事後信用区間

χn2χ^2_nの確率密度関数は

p(y;n)=yn/21ey/22n/2Γ(n/2),y>0. p(y; n) = \frac{y^{n/2 - 1} e^{-y/2}}{2^{n/2} \Gamma(n/2)}, \quad y > 0.

である。y=nνσ2y=\frac{nν}{σ^2}とすると

p(yσ)(nνσ2)n/21enν2σ2(σ2)n/21enν2σ2=σn+2enν2σ2 p(y|σ)∝(\frac{nν}{σ^2})^{n/2 - 1} e^{-\frac{nν}{2σ^2}} \\ ∝(σ^{-2})^{n/2 - 1} e^{-\frac{nν}{2σ^2}} \\ =σ^{-n + 2} e^{-\frac{nν}{2σ^2}}

である。σσの事後分布は

p(σy)p(yσ)p(σ)σn+2enν2σ2σ1=σn+1enν2σ2 p(σ|y)∝p(y|σ)p(σ)∝σ^{-n + 2} e^{-\frac{nν}{2σ^2}}σ^{-1} \\ =σ^{-n + 1} e^{-\frac{nν}{2σ^2}}

σ2σ^2の事後分布は

p(σ2y)p(yσ2)p(σ2)σn+2enν2σ2σ2=σnenν2σ2=(σ2)n/2enν2σ2 p(σ^2|y)∝p(y|σ^2)p(σ^2)∝σ^{-n + 2} e^{-\frac{nν}{2σ^2}}σ^{-2} \\ =σ^{-n} e^{-\frac{nν}{2σ^2}} \\ =(σ^2)^{-n/2} e^{-\frac{nν}{2σ^2}}

である。(a,b)(\sqrt{a},\sqrt{b})σyσ|yの95%最高事後信用区間とするとp(σ=ay)=p(σ=by)p(σ=\sqrt{a}|y)=p(σ=\sqrt{b}|y)となる。

an/2+1/2enν2a=bn/2+1/2enν2b a^{-n/2 + 1/2} e^{-\frac{nν}{2a}}=b^{-n/2 + 1/2} e^{-\frac{nν}{2b}}

(a,b)(\sqrt{a},\sqrt{b})を単純に2乗して変換した(a,b)(a,b)σ2yσ^2|yの95%最高事後信用区間となるか? 95%最高事後信用区間であるならp(σ2=ay)=p(σ2=by)p(σ^2=a|y)=p(σ^2=b|y)となるはずである。

(a)n/2enν2a=(b)n/2enν2b (a)^{-n/2} e^{-\frac{nν}{2a}}=(b)^{-n/2} e^{-\frac{nν}{2b}}

(10)を(11)で割ると

a1/2=b1/2 a^{1/2}=b^{1/2}

となる。したがってa=bである。 しかし、aとbは信用区間の両端なので、これはあり得ない。 よって、(a,b)はσ2yσ^2|yの95%最高事後信用区間ではない。