Probably most common MCMC algorithm and very, very useful. Easily summarized:
z ~ ...
samples = []
for i = 1,...,N:
z' ~ q(z'|z)
x ~ p(x | z); x' ~ p(x'|z')
Q = q(z|z'); Q' = q(z'|z); P = p(x, z); P' = p(x', z')
a = (P' / Q') / (P / Q)
if min(1, a) > uniform():
z = z'
samples.append(z)
@gen function model(n_samples, verbose :: Bool)
b = @trace(beta(2.0, 2.0), :b)
z = @trace(bernoulli(b), :z)
log_sigma = @trace(normal(-1.0, 1.0), :log_sigma)
data = Array{Float64, 1}(undef, n_samples)
for n in 1:n_samples
data[n] = @trace(normal(z, exp(log_sigma)), (:data, n))
end
if verbose
println("Used b = $b, z = $z, and log_sigma = $log_sigma")
end
data
end
DynamicDSLFunction{Any}(Dict{Symbol,Any}(), Dict{Symbol,Any}(), Type[Any, Bool], false, Union{Nothing, Some{Any}}[nothing, nothing], ##model#253, Bool[0, 0], false)
:b
, :z
, and :log_sigma
. We have to be careful to not proposse to any of the (:data, n)
since these are our observable rvs (at least, in this context).mh_step(trace, address) = mh(trace, select(address))[1]
function mh_step(trace, addresses :: AbstractArray)
for address in addresses
trace = mh_step(trace, address)
end
trace
end
function proposal(n_iter, data)
cm = choicemap()
for (n, d) in enumerate(data)
cm[(:data, n)] = d
end
trace, _ = generate(model, (length(data), false), cm)
posterior = Dict(:b => [], :z => [], :log_sigma => [])
addresses = [:b, :z, :log_sigma]
for i in 1:n_iter
trace = mh_step(trace, addresses)
for addr in addresses
push!(posterior[addr], trace[addr])
end
end
posterior
end
proposal (generic function with 1 method)
Random.seed!(2021)
data = model(100, true)
post = proposal(100000, data)
Used b = 0.6105591808638493, z = true, and log_sigma = -1.9209084060204957
Dict{Symbol,Array{Any,1}} with 3 entries: :b => Any[0.488831, 0.807375, 0.75511, 0.75511, 0.418341, 0.898416, 0… :z => Any[false, true, true, true, true, true, true, true, true, true… :log_sigma => Any[-0.902225, -0.902225, -1.29089, -1.91315, -1.91315, -1.9131…
p = plot(p1, p2,)
Random.seed!(2021)
@gen function m1(size :: Int, verbose :: Bool)
loc = @trace(normal(0.0, 1.0), :loc)
log_scale = @trace(normal(0.0, 1.0), :log_scale)
if verbose
println("loc = $loc, log_scale = $log_scale")
end
data = Array{Float64, 1}(undef, size)
for i in 1:size
data[i] = @trace(normal(loc, exp(log_scale)), (:data, i))
end
data
end
DynamicDSLFunction{Any}(Dict{Symbol,Any}(), Dict{Symbol,Any}(), Type[Int64, Bool], false, Union{Nothing, Some{Any}}[nothing, nothing], ##m1#258, Bool[0, 0], false)
# we'll use grad to annotate the fact that we can differentiate w.r.t. parameters in g1
@gen (grad) function g1()
@param loc_mu :: Float64
@param loc_log_sigma :: Float64
@param log_scale_mu :: Float64
@param log_scale_log_sigma :: Float64
@trace(normal(loc_mu, exp(loc_log_sigma)), :loc)
@trace(normal(log_scale_mu, exp(log_scale_log_sigma)), :log_scale)
end
DynamicDSLFunction{Any}(Dict{Symbol,Any}(), Dict{Symbol,Any}(), Type[], false, Union{Nothing, Some{Any}}[], ##g1#259, Bool[], true)
# generate some fake data
size = 500
data = m1(size, true)
data_cm = choicemap()
for (i, pt) in enumerate(data)
data_cm[(:data, i)] = pt
end
# initialize parameters and set optimizer
t_pars = [:loc_mu, :loc_log_sigma, :log_scale_mu, :log_scale_log_sigma] # everything we can grad
for p in t_pars
init_param!(g1, p, 0.0) # more intelligent choices are possible!
end
vi_update = ParamUpdate(GradientDescent(0.0005, 5), g1 => t_pars);
loc = 0.516653904160697, log_scale = -0.08057717729047323
(elbo_estimate, traces, elbo_history) = black_box_vi!(
m1, (size, false),
data_cm,
g1, (),
vi_update,
iters=250
);
plot(p1, p2, p3)