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

Estimation from two independent experiments

ニワトリの脳へのカルシウム流量yを測定する。グループ1はコントロール群で、グループ2は実験群である。

	サンプルサイズ	サンプル平均	サンプル標準偏差
コントロール群	$n_c=32$	$μ_c=1.013$	$σ_c=0.24$
実験群	$n_t=36$	$μ_t=1.173$	$σ_t=0.20$

(a) 平均の事後分布

カルシウム流量yは正規分布に従うとする。その平均と標準偏差は不明である。平均μと標準偏差σの事前分布は一様分布とする。本文p.64より事前分布は

p(μ,σ^2)∝(σ^2)^{-1}

となる。事後分布は

p(μ,σ^2|y)=p(μ,σ^2)p(y|μ,σ^2) \\ =σ^{-n-2}\exp(-\frac{1}{2σ^2}[(n-1)s^2+n(\bar y-μ)^2])

となる。 μの周辺事後分布は

p(μ|y)=\int p(μ,σ^2|y)dσ^2

この計算は本文p.66で行なっている。最終的に

p(µ|y)∝[1 + \frac{n(µ − \bar y)^2}{(n − 1)s^2}]^{−n/2}

となる。これは $t_{n-1}(\bar y,s^2/n)$ 分布である。

µ|y \sim t_{n-1}(\bar y,s^2/n)

今、コントロール群は $N(μ_c,σ_c)$ に、実験群は $N(μ_t,σ_t)$ に従うとしているので、事後分布 $µ_c|y$ 、 $µ_t|y$ はそれぞれ

µ_c|y_c \sim t_{n_c-1}(µ_c,σ_c^2/n_c)=t_{32-1}(1.013,0.24^2/32)

µ_t|y_t \sim t_{n_t-1}(µ_t,σ_t^2/n_t)=t_{36-1}(1.173,0.2^2/36)

である。

using Distributions, StatsPlots
control_dist = 1.013 + TDist(32-1)*0.24/sqrt(32)
exposed_dist = 1.173 + TDist(36-1)*0.2/sqrt(36)
plot(control_dist, label="control", yrotation=90)
plot!(exposed_dist, label="exposed")
savefig(joinpath(@OUTPUT, "TDist.svg"))
nothing

(b) $µ_t-µ_c$ の分布

$µ_t-µ_c$ の分布を調べることで、実験の効果測定を行う。

using Random
Random.seed!(1)
ycontrol=rand(control_dist,5000)
yexposed=rand(exposed_dist,5000)
dif = yexposed .- ycontrol
histogram(dif, label="µt - µc", yrotation=90)
savefig(joinpath(@OUTPUT, "dif.svg"))
quantile(dif,[.025,.975])

2-element Vector{Float64}:
 0.04907497523383009
 0.2686457314755031

図は $µ_t-µ_c$ のヒストグラムであり、95%信用区間は[0.05, 0.27]と推定された。

Estimation from two independent experiments

(a) 平均の事後分布

(b) µt−µcµ_t-µ_cµt​−µc​の分布

(b) $µ_t-µ_c$ の分布