-
Notifications
You must be signed in to change notification settings - Fork 0
/
hmc.jl
62 lines (57 loc) · 1.64 KB
/
hmc.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
doc = """Hamiltonian Monte Carlo (HMC) Sampler.
Usage:
hmc.jl [options] <# of samples>
Options:
-h --help Show this help.
--plot Plot samples.
"""
using Gadfly
using DocOpt
include("util.jl")
# U : potential energy function
# ∇U : gradient of the potential energy function
# θ₀ : initial state
# M : number of samples
# ϵ : step size
# L : number of steps
function hmc(U::Function, ∇U::Function, θ₀::Vector{Float64}, M::Int, ϵ::Float64, L::Int)
d = length(θ₀)
# allocate sampels' holder
samples = Array(typeof(θ₀), M)
# set the current sate to the initail state
θ = θ₀
for m in 1:M
# sample momentum variable
p = randn(d)
H = U(θ) + p ⋅ p / 2
θ̃ = θ
for l in 1:L
p -= ϵ / 2 * ∇U(θ̃) # half step in momentum variable
θ̃ += ϵ * p # full step in location variable
p -= ϵ / 2 * ∇U(θ̃) # half step in momentum variable again
end
H̃ = U(θ̃) + p ⋅ p / 2
if randn() < min(1.0, exp(H - H̃))
# accept the proposal
θ = θ̃
end
samples[m] = θ
print_sample(θ)
end
samples
end
let
args = docopt(doc)
# include functions
include("norm.jl")
M = int(args["<# of samples>"])
for L in [1, 5, 10, 25, 50, 100], ϵ in [0.01, 0.05, 0.1, 0.5]
srand(0)
title = "HMC (ϵ = $ϵ, L = $L)"
println("# $title")
samples = hmc(x -> -ln_normal(x), x -> -∇ln_normal(x), x₀, M, ϵ, L)
if args["--plot"]
plot_samples("$title.svg", samples, title)
end
end
end