Newer
Older
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
using StatsBase, Distributions
using AbstractMCMC: MCMCSerial #, MCMCThreads, MCMCDistributed
using CairoMakie
using Revise
using Pkg; Pkg.activate(".")
using MultilevelChainSampler
function analyse(chains)
t = TruncatedNormal(0, 1, -1, 1)
println(" Target: mean = ", mean(t), ", std = ", std(t))
accept = [ getindex.(c, 1) for c in chains ]
states = [ getindex.(c, 2) for c in chains ]
rejection_rate = 1 .- mean.(accept)
println(" Rejection rate: ", mean(rejection_rate) * 100, " ± ", std(rejection_rate) * 100, " % ")
m = mean.( states )
s = std.( states )
println(" Mean ", mean(m), " ± ", std(m) )
println(" Std ", mean(s), " ± ", std(s) )
println("\n")
fig = Figure(); ax = Axis(fig[1,1])
for s in states[1:1]
hist!(ax, s, normalization=:pdf)
end
x = [-1:.01:1 ... ]; y = pdf.(t, x)
lines!(ax, x, y, label="N(0,1)")
display(fig);
end
begin # Vanilla MetropolisHastings
w = CyclicWalk(-1, 1, .5)
s = MetropolisHastings(w)
f = LogDensity(x -> -x^2/2)
print("Sampling π_∞")
@time chains = sample(f, s, MCMCSerial(), 10000, 4)
analyse(chains)
end
begin # Sample Based Energy Function
w = CyclicWalk(-1, 1, .5)
s = MetropolisHastings(w)
n = 100
z = [ rand(2) .* 2 .- 1 for i=1:n ]
f = SampledLogDensity(z, (x,z) -> - ( sum(z.^2) < x^2 )/2 )
print("Sampling π_$n")
@time chains = sample(f, s, MCMCSerial(), 10000, 4)
analyse(chains)
end