Exercise 5.1 Solution Example - Hoff, A First Course in Bayesian Statistical Methods
標準ベイズ統計学 演習問題 5.1 解答例
a)
answer
school1 = readdlm("./Exercises/school1.dat")
school2 = readdlm("./Exercises/school2.dat")
school3 = readdlm("./Exercises/school3.dat")
# set prior parameters
μ₀ = 5
σ₀² = 4
κ₀ = 1
ν₀ = 2
# %%
n1 = length(school1)
n2 = length(school2)
n3 = length(school3)
ȳ₁ = mean(school1)
ȳ₂ = mean(school2)
ȳ₃ = mean(school3)
s₁² = var(school1)
s₂² = var(school2)
s₃² = var(school3)
function calc_posterior_param(n, ȳ, s², μ₀, σ₀², κ₀, ν₀)
κn = κ₀ + n
νn = ν₀ + n
μn = (κ₀ * μ₀ + n * ȳ) / κn
σ²n = (ν₀ * σ₀² + (n-1) * s² + (κ₀ * n * (ȳ - μ₀)^2)/ κn ) / νn
return κn, νn, μn, σ²n
end
function calc_mc(n, ȳ, s², μ₀, σ₀², κ₀, ν₀, n_mc)
κn, νn, μn, σ²n = calc_posterior_param(n, ȳ, s², μ₀, σ₀², κ₀, ν₀)
dist_inverse_σ² = Gamma(νn/2, 2/(νn * σ²n))
σ²_mc = rand(dist_inverse_σ², n_mc) .^ -1
calc_θ_mc = σ² -> rand(Normal(μn, sqrt(σ²/κn)))
θ_mc = calc_θ_mc.(σ²_mc)
return θ_mc, sqrt.(σ²_mc)
end
function write_quantile_ci(θ_mc, σ_mc, ci)
lower = (1 - ci) / 2
upper = 1 - lower
θ_quantile = round.(quantile(θ_mc, [lower, upper]), digits=3)
σ_quantile = round.(quantile(σ_mc, [lower, upper]), digits=3)
res = ["$(θ_quantile[1]) < θ < $(θ_quantile[2])"
,"$(σ_quantile[1]) < σ < $(σ_quantile[2])"]
return res
end
# School 1
θ₁_mc, σ₁_mc = calc_mc(n1, ȳ₁, s₁², μ₀, σ₀², κ₀, ν₀, 5000);
write_quantile_ci(θ₁_mc, σ₁_mc, 0.95)
2-element Vector{String}:
"7.746 < θ < 10.814"
"2.985 < σ < 5.181"
# School 2
θ₂_mc, σ₂_mc = calc_mc(n2, ȳ₂, s₂², μ₀, σ₀², κ₀, ν₀, 5000);
write_quantile_ci(θ₂_mc, σ₂_mc, 0.95)
2-element Vector{String}:
"5.24 < θ < 8.753"
"3.356 < σ < 5.91"
# School 3
θ₃_mc, σ₃_mc = calc_mc(n3, ȳ₃, s₃², μ₀, σ₀², κ₀, ν₀, 5000);
write_quantile_ci(θ₃_mc, σ₃_mc, 0.95)
2-element Vector{String}:
"6.129 < θ < 9.421"
"2.818 < σ < 5.114"
b)
answer
θ_mcs = [θ₁_mc, θ₂_mc, θ₃_mc]
perms(l) = isempty(l) ? [l] : [[x; y] for x in l for y in perms(setdiff(l, x))]
for perm in perms([1,2,3])
i, j, k = perm
p = mean(θ_mcs[i] .< θ_mcs[j] .< θ_mcs[k])
println("p(θ$i < θ$j < θ$k) = $p")
end
p(θ1 < θ2 < θ3) = 0.0084 p(θ1 < θ3 < θ2) = 0.005 p(θ2 < θ1 < θ3) = 0.0904 p(θ2 < θ3 < θ1) = 0.6506 p(θ3 < θ1 < θ2) = 0.02 p(θ3 < θ2 < θ1) = 0.2256
c)
answer
function sample_posterior_predictive(θ_mc, σ_mc)
ỹ_mc = rand.(Normal.(θ_mc, σ_mc))
return ỹ_mc
end
ỹ₁_mc = sample_posterior_predictive(θ₁_mc, σ₁_mc)
ỹ₂_mc = sample_posterior_predictive(θ₂_mc, σ₂_mc)
ỹ₃_mc = sample_posterior_predictive(θ₃_mc, σ₃_mc)
ỹ_mcs = [ỹ₁_mc, ỹ₂_mc, ỹ₃_mc]
for perm in perms([1,2,3])
i, j, k = perm
p = mean(ỹ_mcs[i] .< ỹ_mcs[j] .< ỹ_mcs[k])
println("p(ỹ$i < ỹ$j < ỹ$k) = $p")
end
p(ỹ1 < ỹ2 < ỹ3) = 0.109 p(ỹ1 < ỹ3 < ỹ2) = 0.1008 p(ỹ2 < ỹ1 < ỹ3) = 0.1882 p(ỹ2 < ỹ3 < ỹ1) = 0.2634 p(ỹ3 < ỹ1 < ỹ2) = 0.1498 p(ỹ3 < ỹ2 < ỹ1) = 0.1888
d)
answer
Using the results from the previous calculations,
\begin{align*} Pr(\theta_1 > \theta_2, \theta_1 > \theta_3 | \text{data}) &= Pr(\theta_1 > \theta_2 > \theta_3 | \text{data}) + Pr(\theta_1 > \theta_3 > \theta_2 | \text{data}) \\ &= 0.6506 + 0.2256 = 0.8762 \\ \\ Pr(\tilde{Y}_1 > \tilde{Y}_2, \tilde{Y}_1 > \tilde{Y}_3 | \text{data}) &= Pr(\tilde{Y}_1 > \tilde{Y}_2 > \tilde{Y}_3 | \text{data}) + Pr(\tilde{Y}_1 > \tilde{Y}_3 > \tilde{Y}_2 | \text{data}) \\ &= 0.2634 + 0.1888 = 0.4522 \end{align*}
