Juliaの世界

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

Posterior distribution as a compromise between prior information and data

コインをn回振る。yは表が出た回数。表が出る確率はθ。

(a) yの事前予測分布

yの事前予測分布Pr(y=k)を計算せよ。

θは一様分布Uniform(0,1)に従うとする。yの事前予測分布は

Pr(y=k)=01Pr(y=kθ)Pr(θ)dθ Pr(y=k)=\int^1_0 Pr(y=k|θ)Pr(θ)dθ

で計算できる。Pr(θ)=1なので

Pr(y=k)=01Pr(y=kθ)dθ Pr(y=k)=\int^1_0 Pr(y=k|θ)dθ

yは二項分布Binomial(n,θ)に従うので、そのPDFは

Pr(y=kθ)=(nk)θk(1θ)nk Pr(y=k|θ)=\binom{n}{k}θ^k(1-θ)^{n-k}

である。従って

Pr(y=k)=01(nk)θk(1θ)nkdθ Pr(y=k)=\int^1_0 \binom{n}{k}θ^k(1-θ)^{n-k}dθ

となる。ここで、ベータ関数は

B(z1,z2)=01tz11(1t)z21 dt \mathrm {B} (z_{1},z_{2})=\int _{0}^{1}t^{z_{1}-1}(1-t)^{z_{2}-1}\,dt

であるので、

01θ(k+1)1(1θ)(nk+1)1dθ=Γ(k+1) Γ(nk+1)Γ(k+1+nk+1)=Γ(k+1) Γ(nk+1)Γ(n+2) \int _{0}^{1}θ^{(k+1)-1}(1-θ)^{(n-k+1)-1}dθ \\ =\frac {\Gamma (k+1)\,\Gamma (n-k+1)}{\Gamma (k+1+n-k+1)} \\ =\frac {\Gamma (k+1)\,\Gamma (n-k+1)}{\Gamma (n+2)}

である。ベータ関数は

B(a,b)=Γ(a) Γ(b)Γ(a+b) B(a,b)=\frac {\Gamma (a)\,\Gamma (b)}{\Gamma (a+b)}

とガンマ関数を用いて表現できる。ガンマ関数は整数の範囲で

Γ(k)=(k1)! \Gamma (k)=(k-1)!

となる。全部まとめて、

Pr(y=k)=n!k!(nk)!k!(nk)!(n+1)!=1n+1 Pr(y=k)=\frac{n!}{k!(n-k)!}\frac{k!(n-k)!}{(n+1)!} \\ =\frac{1}{n+1}

である。 コインの表が出る確率θを0から1で一様と考えているとき、表が出る回数yの事前予測分布はコインを振る回数nの範囲内で一様分布している。

(b) θの事後分布の平均

θはBeta(α,β)に従うとする。今、n回コインを振り、y回表が出たことを観測した。 θの事後分布の平均α+yα+β+n\frac{α+y}{α+β+n}は常に事前分布の平均αα+β\frac{α}{α+β}と相対頻度yn\frac{y}{n}の間に位置することを確認せよ。

zがxとyの間に位置するという関係は1次元の凸結合を考えればよい。 z=ax+byでa+b=1の時、zはxとyの線分上にある。 事後分布の平均を以下の関係式で表現する。

α+yα+β+n=tαα+β+(1t)yn \frac{α+y}{α+β+n}=t\frac{α}{α+β}+(1-t)\frac{y}{n}

0≤t≤1であれば、その2つの間にあるということになる。tに関して解くと

t=α+βα+β+n t=\frac{α+β}{α+β+n}

となる。よってtは常に0<t<1であるので、 θの事後分布の平均α+yα+β+n\frac{α+y}{α+β+n}は常に事前分布の平均αα+β\frac{α}{α+β}と相対頻度yn\frac{y}{n}の間に位置する。

(c) θの事後分布の分散

θは一様分布に従うとする。θの事後分布の分散は常に事前分布の分散より小さいことを示せ。

Beta(1,1)はUniform(0,1)と同等である。 Beta(α,β)の分散はαβ(α+β)2(α+β+1)\frac{αβ}{(α+β)^2(α+β+1)}である。 よって、θの事前分布の分散は122(1+1+1)=112\frac{1}{2^2(1+1+1)}=\frac{1}{12}となる。 Beta分布は二項分布の共役事前分布なのでθの事後分布はBeta(1+y,1+n-y)となる。 その分散は

αβ(α+β)2(α+β+1)=(1+y)(1+ny)(1+y+1+ny)2(1+y+1+ny+1)=(1+y)(1+ny)(2+n)2(3+n) \frac{αβ}{(α+β)^2(α+β+1)}=\frac{(1+y)(1+n-y)}{(1+y+1+n-y)^2(1+y+1+n-y+1)} \\ =\frac{(1+y)(1+n-y)}{(2+n)^2(3+n)}

相加相乗平均の不等式xyx+y2\sqrt{xy}≤\frac{x+y}{2}xy14(x+y)2xy≤\frac{1}{4}(x+y)^2を使う。

14(1+y2+n+1+ny2+n)13+n=1413+n<112 ≤\frac{1}{4}(\frac{1+y}{2+n}+\frac{1+n-y}{2+n})\frac{1}{3+n} \\ =\frac{1}{4}\frac{1}{3+n} \\ <\frac{1}{12}

従って、θの事後分布の分散は事前分布の分散1/12より常に小さい。

(d) θの事後分布の分散の確認

(c)は一様分布(Beta(1,1))の時のみ分散が小さくなることを確認したが、他の場合はどうなるか確認する。

using Distributions
var.([Beta(1,1), Beta(2,2), Beta(3,5), Beta(5,3)])
4-element Vector{Float64}:
 0.08333333333333333
 0.05
 0.026041666666666668
 0.026041666666666668

この結果から言えることはBeta分布の分散はα、βが大きくなれば常に小さくなるということである。 事前分布Beta(a,b)に対し、事後分布はBeta(a+y,b+n-y)であるので、事後分布のα、βは常に大きくなる。 そのため、事後分布の分散は事前分布の分散より常に小さい。