Exercise 10.2 Solution Example - Hoff, A First Course in Bayesian Statistical Methods
標準ベイズ統計学 演習問題 10.2 解答例
a)
Answer
b)
Answer
Let \(p\) be the probability of nesting success. Then,
\begin{align*} \alpha + \beta x_i &= \text{logit } p \\ &= \log \frac{p}{1-p}. \\ \end{align*}Thus, if we have no prior information about the probability of nesting success, we can set the prior expectation of \(\alpha + \beta x_i\) to 0, \(p = 0.5\). Also, we have no prior information about the relationship between wingspan and nesting success, so we can set the prior expectation of \(\beta\) to 0. Consequently, we can set the prior expectation of \(\alpha\) to 0.
Next, we consider the prior variance of \(\alpha\) and \(\beta\). Suppose temporarily that \(\beta = 0\). If the prior variance of \(\alpha\) is set to \(2500\), \(\alpha\) ranges between \(-100\) and \(100\) for 95% probability, and the odds ratio should be between \(\exp(-100) \approx 3.7 \times 10^{-44}\) and \(\exp(100) \approx 2.7 \times 10^{43}\). This prior can be seen quite uninformative. Similarly, considering \(x_i\) ranges between 10 and 15, we can set the prior variance of \(\beta\) to 25. Thus, we can set the prior distribution of \(\alpha\) and \(\beta\) as follows:
\begin{align*} \alpha &\sim \mathcal{N}(0, 2500) \\ \beta &\sim \mathcal{N}(0, 25) \\ \end{align*}c)
Answer
using BayesPlots
using DelimitedFiles
using Distributions
using LaTeXStrings
using LinearAlgebra
using MCMCChains
using Plots
using Plots.Measures
using Random
using StatsBase
using StatsModels
using StatsPlots
using GLM
data = readdlm("../../Exercises/msparrownest.dat")
y = data[:, 1]
x = data[:, 2]
# add intercept
x = hcat(ones(length(x)), x)
function MetropolisLogit(S, y, x, β₀, Σ₀, var_prop; seed=123)
# S: number of samples
# y: response vector
# x: predictor vector
# β₀: prior mean of β
# Σ₀: prior covariance of β
# var_prop: proposal variance
function logistic(ψ)
return exp(ψ) / (1 + exp(ψ))
end
Random.seed!(seed)
q = length(β₀)
BETA = Matrix{Float64}(undef, S, q)
acs = 0 # acceptance rate
# Initialize
glmfit = glm(x, y, Binomial(), LogitLink())
β = coef(glmfit)
# Metropolis algorithm
for s in 1:S
p = logistic.(x*β)
# Draw β*
β_star = rand(MvNormal(β, var_prop))
p_star = logistic.(x*β_star)
log_r = sum(logpdf.(Bernoulli.(p_star), y)) +
logpdf(MvNormal(β_star, Σ₀), β) -
sum(logpdf.(Bernoulli.(p), y)) -
logpdf(MvNormal(β, Σ₀), β_star)
if log(rand()) < log_r
β = β_star
acs += 1
end
BETA[s, :] = β
end
return BETA, acs/S
end
# Prior
β₀ = [0, 0]
Σ₀ = [2500 0; 0 25]
S = 10000
# search best proposal variance
for i in 1:20
# Proposal variance
var_prop = i*inv(x'x)
BETA, acs = MetropolisLogit(S, y, x, β₀, Σ₀, var_prop)
println("i = ", i)
println("acceptance rate = ", acs)
ess_α = MCMCDiagnosticTools.ess(BETA[:,1])
ess_β = MCMCDiagnosticTools.ess(BETA[:,2])
println("ess_α = ", ess_α)
println("ess_β = ", ess_β)
end
var_prop = 17*inv(x'x)
BETA, acs = MetropolisLogit(S, y, x, β₀, Σ₀, var_prop)
ess(MCMCChains.Chains(BETA))
ESS
parameters ess ess_per_sec
Symbol Float64 Missing
param_1 1261.2209 missing
param_2 1248.0384 missing
d)
Answer
e)
Answer
x_grid = 10:0.1:15
x_grid = hcat(ones(length(x_grid)), x_grid)
# %%
f_post = logistic.(x_grid*BETA[burnin:end, :]')
# %%
f_post_mean = mean(f_post, dims=2)
f_post_lower = map(x -> quantile(x, 0.025), eachrow(f_post))
f_post_upper = map(x -> quantile(x, 0.975), eachrow(f_post))
