Exercise 5.1 Solution Example - Hoff, A First Course in Bayesian Statistical Methods
標準ベイズ統計学演習問題 5.1 解答例

a)
- answer
b)
- answer
c)
- answer
d)
- answer

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)

ȳ₁ = mean(school1)
ȳ₂ = mean(school2)
ȳ₃ = mean(school3)

s₁² = var(school1)
s₂² = var(school2)
s₃² = var(school3)

function calc_posterior_param(n, ȳ, s², μ₀, σ₀², κ₀, ν₀)
    κn = κ₀ + n
    νn = ν₀ + n
    μn =  (κ₀ * μ₀ + n * ȳ) / κn
    σ²n = (ν₀ * σ₀² + (n-1) * s² + (κ₀ * n * (ȳ - μ₀)^2)/ κn ) / νn

    return κn, νn, μn, σ²n
end

function calc_mc(n, ȳ, s², μ₀, σ₀², κ₀, ν₀, n_mc)
    κn, νn, μn, σ²n = calc_posterior_param(n, ȳ, 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, ȳ₁, 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, ȳ₂, 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, ȳ₃, 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)
    ỹ_mc = rand.(Normal.(θ_mc, σ_mc))
    return ỹ_mc
end

ỹ₁_mc = sample_posterior_predictive(θ₁_mc, σ₁_mc)
ỹ₂_mc = sample_posterior_predictive(θ₂_mc, σ₂_mc)
ỹ₃_mc = sample_posterior_predictive(θ₃_mc, σ₃_mc)

ỹ_mcs = [ỹ₁_mc, ỹ₂_mc, ỹ₃_mc]

for perm in perms([1,2,3])
    i, j, k = perm
    p = mean(ỹ_mcs[i] .< ỹ_mcs[j] .< ỹ_mcs[k])
    println("p(ỹ$i < ỹ$j < ỹ$k) = $p")
end

p(ỹ1 < ỹ2 < ỹ3) = 0.109
p(ỹ1 < ỹ3 < ỹ2) = 0.1008
p(ỹ2 < ỹ1 < ỹ3) = 0.1882
p(ỹ2 < ỹ3 < ỹ1) = 0.2634
p(ỹ3 < ỹ1 < ỹ2) = 0.1498
p(ỹ3 < ỹ2 < ỹ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*}

4.8

5.2

Exercise 5.1 Solution Example - Hoff, A First Course in Bayesian Statistical Methods 標準ベイズ統計学 演習問題 5.1 解答例

Table of Contents

a)

answer

b)

answer

c)

answer

d)

answer

Exercise 5.1 Solution Example - Hoff, A First Course in Bayesian Statistical Methods
標準ベイズ統計学演習問題 5.1 解答例