mschauer
commented
3 years ago Example:
using ZigZagBoomerang, Distributions, ForwardDiff, LinearAlgebra, SparseArrays, StructArrays
using Random
d = 2
ℓ(x) = logpdf(MvNormal([1.0, 2.0], Diagonal([0.1, 10.0]).^2), x)
function negpartiali(f,d)
id = collect(I(d))
ith(i) = @inbounds @view id[:,i]
function (x,i)
sa = StructArray{ForwardDiff.Dual{}}((x, ith(i)))
δ = f(sa).partials[]
return -δ
end
end
∇ϕi2 = negpartiali(ℓ, 2)
∇ϕi_randomized(x, i) = ∇ϕi2(x, i) + 0.5*(randexp() - 1.0)
t0 = 0.0
x0 = rand(d)
θ0 = 1.0ones(d)
c= 500.0
T = 5000.0
Z = ZigZag(sparse(I(d)), 0*x0, [1.0, 1.0], 5.0, 0., 1.0)
trace, final, (acc, num) = spdmp(∇ϕi_randomized, t0, x0, θ0, T, c*ones(d), Z; adapt=true, adaptscale=true, progress=true)
ts, xs = ZigZagBoomerang.sep(discretize(trace, 0.1))
a = [xs[i][1] for i in 1:length(ts)]
b = [xs[i][2] for i in 1:length(ts)]
mean(trace)
@show Z.σ, (/(Z.σ...))^2
@show acc/T
r = eachindex(ts)[1:end]
using GLMakie
fig = Figure()
ax = fig[1,1] = Axis(fig)
lines!(ax, first.(xs), last.(xs))
ax = fig[1,2] = Axis(fig)
lines!(ax, ts[r], first.(xs)[r])
ax = fig[1,3] = Axis(fig)
lines!(ax, ts[r],(last.(xs))[r])
fig
Using
primal/dual averaging of NUTS