Independent MH
Gamma distribution
Robert and Casella (2013) presents the following algorithm:

If there exists a constant such that
the algorithm produces a uniformly ergodic chain (Theorem), and the expected acceptance probability associated with the algorithm is at least when the chain is stationary, and in that sense, the IMH is more efficient than the Accept-Reject algorithm.
Let's illustrate this algorithm with . We have introduced how to sample from Gamma distribution via Accept-Reject algorithm in Special Distributions, and it is straightforward to get the Gamma Metropolis-Hastings based on the ratio of ,

And we can implement this algorithm with the Julia code:
function mh_gamma(T = 100, alpha = 1.5)
a = Int(floor(alpha))
b = a/alpha
x = ones(T+1) # initial value: 1
for t = 1:T
yt = rgamma_int(a, b)
rt = (yt / x[t] * exp((x[t] - yt) / alpha))^(alpha-a)
if rt >= 1
x[t+1] = yt
else
u = rand()
if u < rt
x[t+1] = yt
else
x[t+1] = x[t]
end
end
end
return(x)
end
To sample , and estimate .
# comparison with accept-reject
res = mh_gamma(5000, 2.43)[2:end]
est = cumsum(res.^2) ./ collect(1:5000)
res2 = ones(5000)
for i = 1:5000
res2[i] = rgamma(2.43, 1)
end
est2 = cumsum(res2.^2) ./ collect(1:5000)
using Plots
plot(est, label="Independent MH")
plot!(est2, label="Accept-Reject")
hline!([8.33], label="True value")

Logistic Regression
We observe according to the model
The likelihood is
and let and put a flat prior on , i.e.,
Note that
where
is the Euler's Constant.
Choose the data-dependent value that makes , where is the MLE of , so .
We can use MLE to estimate the coefficient in the logistical model (see my post for the derivation), and the following Julia code can help us fit the model quickly.
using DataFrames, GLM, Plots
temp = [53, 57, 58, 63, 66, 67, 67, 67, 68, 69, 70, 70, 70, 70, 72, 73, 75, 75, 76, 76, 78, 79, 81]
failure = [1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0]
data = DataFrame(temp = temp, failure = failure)
logit_fit = glm(@formula(failure ~ temp), data, Binomial(), LogitLink())
plot(temp, predict(logit_fit), legend = false, xlabel = "Temperature", ylab = "Probability")
scatter!(temp, predict(logit_fit))

The estimates of the parameters are and .
Wrote the following Julia code to implement the independent MH algorithm,
## metropolis-hastings
using Distributions
γ = 0.57721
function ll(α::Float64, β::Float64)
a = exp.(α .+ β*temp)
return prod( (a ./ (1 .+ a) ).^failure .* (1 ./ (1 .+ a)).^(1 .- failure) )
end
function mh_logit(T::Int, α_hat::Float64, β_hat::Float64, σ_hat::Float64)
φ = Normal(β_hat, σ_hat)
π = Exponential(exp(α_hat+γ))
Α = ones(T)
Β = ones(T)
for t = 1:T-1
α = log(rand(π))
β = rand(φ)
r = ( ll(α, β) / ll(Α[t], Β[t]) ) * ( pdf(φ, β) / pdf(φ, Β[t]) )
if rand() < r
Α[t+1] = α
Β[t+1] = β
else
Α[t+1] = Α[t]
Β[t+1] = Β[t]
end
end
return Α, Β
end
Α, Β = mh_logit(10000, 15.04, -0.233, 0.108)
p1 = plot(Α, legend = false, xlab = "Intercept")
hline!([15.04])
p2 = plot(Β, legend = false, xlab = "Slope")
hline!([-0.233])
plot(p1, p2, layout = (1,2))

Saddlepoint tail area approximation
If is the cumulant generating function, solving the saddlepoint equation yields the saddlepoint. For noncentral chi squared random variable , the moment generating function is
where is the number of degrees of freedom and is the noncentrality parameter, and its saddlepoint is
The saddlepoint can be used to approximate the tail area of a distribution. We have the approximation
where is the sample from the saddlepoint distribution. Using a Taylor series approximation,
so the instrumental density can be chosen as , where
I implemented the independent MH algorithm via the following code to produce random variables from the saddlepoint distribution.
using Distributions
p = 6
λ = 9
function K(t)
return 2λ*t / (1-2t) - p / 2 * log(1-2t)
end
function D1K(t)
return ( 4λ*t ) / (1-2t)^2 + ( 2λ + p ) / ( 1 - 2t )
end
function D2K(t)
return 2(p*(1-2t) + 4λ) / (1-2t)^3
end
function logf(n, t)
return n * (K(t) - t*D1K(t)) + 0.5log(D2K(t))
end
function logg(n, t)
return -1n * D2K(0) * t^2/2
end
function mh_saddle(T::Int = 10000; n::Int = 1)
Z = zeros(T)
g = Normal(0, 1 / sqrt(n*D2K(0)))
for t = 1:T-1
z = rand(g)
logr = logf(n, z) - logf(n, Z[t]) + logg(n, Z[t]) - logg(n, z)
if log(rand()) < logr
Z[t+1] = z
else
Z[t+1] = Z[t]
end
end
return Z
end
Z = mh_saddle(n = 1)
function tau(x)
return ( -1p + 2x - sqrt(p^2 + 8λ * x) ) / (4x)
end
println(sum(Z .> tau(36.225)) / 10000)
println(sum(Z .> tau(40.542)) / 10000)
println(sum(Z .> tau(49.333)) / 10000)
I can successfully reproduce the results.
If ,
Interval
1 - pchisq(x, 6, 9*2)
IMH
0.1000076
0.1037
0.05000061
0.0516
0.01000057
0.0114
For , note that if , then , then we can use the following R code to figure out the cutpoints.
qchisq(0.90, 600, 1800)/100
qchisq(0.95, 600, 1800)/100
qchisq(0.99, 600, 1800)/100
Then use these cutpoints to estimate the probability by using samples produced from independent MH algorithm.
Interval
1 - pchisq(x*100, 6*100, 9*2*100)
IMH
0.10
0.1045
0.05
0.0516
0.01
0.0093
Last updated