ARMS MH

Search…
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!
Random Walk MH
PBMCMC
Last modified 2yr ago
Copy link