ARMS MH
ARMS stands for Adaptive Rejection Metropolis Sampling, and it is the generalization of ARS algorithm.
We can implement this algorithm with the following Julia program:
1
include("../GenRV/ars.jl")
2
using Main.corears
3
# ARMS
4
function arms(T, yfixed::Array)
5
x = ones(T+1)
6
for t = 1:T
7
# generate Y
8
while true
9
global y
10
y = gplus_sample(yfixed)
11
u = rand()
12
u <= exp(h(y)-hplus(y, yfixed)) && break
13
end
14
# accept or not
15
v = rand()
16
r = exp(h(y))*phi(x[t], yfixed)/(exp(h(x[t]))*phi(y, yfixed))
17
println(r)
18
if r >= 1
19
x[t+1] = y
20
else
21
if v <= r
22
x[t+1] = y
23
else
24
x[t+1] = x[t]
25
end
26
end
27
end
28
return(x)
29
end
Copied!
Now let's apply ARMS to poisson logistic model:
1
# data
2
x = rand(10)
3
y = rand(10)
4
# parameters
5
tau = 1
6
7
# function of h
8
function h(b)
9
return(b * sum(x .* y) - sum(log.(1 .+ exp.(b*x))) - sum(1 .- 1 ./ (1 .+ exp.(b*x))) - b^2/(2*tau^2))
10
end
11
12
# derivative of h
13
function dh(b)
14
return(sum(x .* y) - sum(x .* (1 .- 1 ./ (1 .+ exp.(b*x)) )) + sum(x .* exp.(b*x) ./ (1 .+ exp.(b*x)).^2 ) - b/tau^2)
15
end
16
17
function phi(x, yfixed::Array)
18
if exp(h(x)) <= exp(hplus(x, yfixed))
19
return(exp(h(x)))
20
else
21
return(exp(hplus(x, yfixed)))
22
end
23
end
24
## example
25
x = arms(10000, [-1.5, -0.7, 0.7, 1.5])
Copied!
Last modified 2yr ago
Copy link