Exercise 8.3 Solution Example - Hoff, A First Course in Bayesian Statistical Methods
標準ベイズ統計学 演習問題 8.3 解答例
a)
Answer
using Distributions using Random using LaTeXStrings using Plots, StatsPlots using Plots.Measures using CSV, DataFrames using DelimitedFiles using StatsBase using SpecialFunctions Random.seed!(1234) # %% # read data Y = Matrix{Float64}(undef, 0, 2) for i in 1:8 school = readdlm("../../Exercises/school$i.dat") Y = vcat(Y, hcat(repeat([i], size(school, 1)), school)) end # %% # set prior parameters μ₀ = 7 γ₀² = 5 τ₀² = 10 η₀ = 2 σ₀² = 10 ν₀ = 2 # %% # starting values m = length(unique(Y[:,1])) n = Vector{Float64}(undef, m) ȳ = Vector{Float64}(undef, m) sv = Vector{Float64}(undef, m) for j in 1:m ȳ[j] = mean(Y[Y[:,1] .== j, 2]) sv[j] = var(Y[Y[:,1] .== j, 2]) n[j] = length(Y[Y[:,1] .== j, 2]) end θ = copy(ȳ) σ² = mean(sv) μ = mean(θ) τ² = var(θ) # setup MCMC Random.seed!(1) S = 5000 THETA = Matrix{Float64}(undef, S, m) MST = Matrix{Float64}(undef, S, 3) # MCMC algorithm for s in 1:S # sample new values of the θs for j in 1:m v_θ = 1 / (n[j] / σ² + 1 / τ²) E_θ = v_θ * (n[j]*ȳ[j]/σ² + μ/τ²) θ[j] = rand(Normal(E_θ, sqrt(v_θ))) end # sample new value of σ² νₙ = ν₀ + sum(n) ss = ν₀ * σ₀² for j in 1:m ss += sum((Y[Y[:,1] .== j, 2] .- θ[j]).^2) end σ² = 1 / rand(Gamma(νₙ / 2, 2/ss)) # sample a new value of μ v_μ = 1 / (m/τ² + 1/γ₀²) E_μ = v_μ * (m * mean(θ) / τ² + μ₀ / γ₀²) μ = rand(Normal(E_μ, sqrt(v_μ))) # sample a new value of τ² etam = η₀ + m ss = η₀ * τ₀² + sum((θ .- μ).^2) τ² = 1 / rand(Gamma(etam/2, 2/ss)) # store results THETA[s, :] = θ MST[s, :] = [μ, σ², τ²] end
Figure 1: stationarity plot
using MCMCDiagnosticTools using MCMCChains ess(MCMCChains.Chains(MST))
ESS
parameters ess ess_per_sec
Symbol Float64 Missing
param_1 4046.5607 missing
param_2 4697.5784 missing
param_3 3299.9524 missing
b)
Answer
# posterior mean and 95% confidence region for σ² julia> mean(MST[:,2]) 14.413317696499103 julia> quantile(MST[:,2], [0.025, 0.975]) 2-element Vector{Float64}: 11.733606090321313 17.92710379952854
# posterior mean and 95% confidence region for μ julia> mean(MST[:,1]) 7.587311747784735 julia> quantile(MST[:,1], [0.025, 0.975]) 2-element Vector{Float64}: 5.910059688648781 9.160791496009686
julia> mean(MST[:,3]) 5.478604048700004 julia> quantile(MST[:,3], [0.025, 0.975]) 2-element Vector{Float64}: 1.9240022529423977 14.101516261365546
Figure 2: Comparison of posterior and prior densities
The following things were learned from the data:
- \(\sigma^2\) is around 14.4, which is higher than the prior mean 10.0, and the variance is smaller than the prior variance 100.
- \(\mu\) is around 7.6, which is higher than the prior mean 7.0, and the variance is smaller than the prior variance 10.0.
- \(\tau^2\) is around 5.5, which is smaller than the prior mean 10.0, and the variance is smaller than the prior variance 100.
c)
Answer
σ²_pri = rand(Gamma(ν₀/2, 2σ₀²/ν₀), S) τ²_pri = rand(Gamma(η₀/2, 2τ₀²/η₀), S) R_pri = τ²_pri ./ (σ²_pri .+ τ²_pri) R_pos = MST[:,3] ./(MST[:,2] .+ MST[:,3])
Figure 3: Prior and posterior densities of \(R\)
The plot above shows that about 25% of the variance is due to between-school variation, and the rest is due to within-school variation.
d)
Answer
# posterior probability that θ₇ is smaller than θ₆ mean(THETA[:,7] .< THETA[:,6])
0.5222
# posterior probability that θ₇ is the smallest of all the θ’s mean(THETA[:,7] .≤ minimum(THETA, dims=2))
0.325
e)
Answer
θ̄ = mean(THETA, dims=1) |> vec fig_e = plot( ȳ, θ̄, seriestype=:scatter, xlabel=L"\bar{y}_i", ylabel=L"\theta_i", legend=:none, xlims=(6,11), ylims=(6,11), size=(400, 400), margin=5mm, alpha=0.5 ) fig_e = plot!(fig_e, x->x, color=:red, linewidth=1,alpha=0.5) fig_e
Figure 4: Sample averages and posterior expectations
For the schools with large sample averages, the posterior expectations tend to be smaller than the sample averages. For the schools with small sample averages, the posterior expectations tend to be larger than the sample averages.
# sample mean of all observations mean(Y[:,2])
7.6912777777777785
# posterior mean of μ mean(MST[:,1])
7.587311747784735
