Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix HMC in Gibbs issue with dynamic models #1217

Merged
merged 6 commits into from
Apr 13, 2020
Merged

Fix HMC in Gibbs issue with dynamic models #1217

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
3c2692b
fix #1095
mohamed82008 Apr 13, 2020
ccbe2d3
add test
mohamed82008 Apr 13, 2020
bf18a62
fix bug
mohamed82008 Apr 13, 2020
ab21f01
test PG and MH
mohamed82008 Apr 13, 2020
c8adde2
fix tests
mohamed82008 Apr 13, 2020
c94e928
avoid duplicate line
mohamed82008 Apr 13, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 9 additions & 6 deletions src/inference/hmc.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -337,23 +337,26 @@ function AbstractMCMC.step!(
spl.state.eval_num = 0

Turing.DEBUG && @debug "current ϵ: $ϵ"
# Gibbs component specified cares

# When a Gibbs component
if spl.selector.tag != :default
# Transform the space
Turing.DEBUG && @debug "X-> R..."
link!(spl.state.vi, spl)
runmodel!(model, spl.state.vi, spl)
end
# Get position and log density before transition
θ_old, log_density_old = spl.state.vi[spl], getlogp(spl.state.vi)
if spl.selector.tag != :default
# Update Hamiltonian
metric = gen_metric(length(spl.state.vi[spl]), spl)
metric = gen_metric(length(θ_old), spl)
∂logπ∂θ = gen_∂logπ∂θ(spl.state.vi, spl, model)
logπ = gen_logπ(spl.state.vi, spl, model)
spl.state.h = AHMC.Hamiltonian(metric, logπ, ∂logπ∂θ)
resize!(spl.state.z.θ, length(θ_old))
Copy link
Member

@yebai yebai Apr 13, 2020 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is worth noting that the vanilla HMC with identity mass matrix and pre-specified leapfrog step-size might be the only sensible HMC sampler for dynamic models for now, since the adaption schedules for step-size and mass matrix are not guaranteed to work well when the static dimensionality condition no longer holds.

@trappmartin @mohamed82008 @xukai92

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@trappmartin Maybe add a warning for this in the Infinite Mixture Model example?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good point. I'll open a PR for this.

spl.state.z.θ .= θ_old
end

# Get position and log density before transition
θ_old, log_density_old = spl.state.vi[spl], getlogp(spl.state.vi)

# Transition
t = AHMC.step(rng, spl.state.h, spl.state.traj, spl.state.z)
# Update z in state
Expand Down
26 changes: 26 additions & 0 deletions test/inference/gibbs.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@ using Random, Turing, Test
import AbstractMCMC
import MCMCChains
import Turing.Inference
using Turing.RandomMeasures

dir = splitdir(splitdir(pathof(Turing))[1])[1]
include(dir*"/test/test_utils/AllUtils.jl")
Expand Down Expand Up @@ -111,4 +112,29 @@ include(dir*"/test/test_utils/AllUtils.jl")
alg = Gibbs(MH(:s), HMC(0.2, 4, :m))
sample(model, alg, 100; callback = callback)
end

@turing_testset "dynamic model" begin
@model imm(y, alpha, ::Type{M}=Vector{Float64}) where {M} = begin
N = length(y)
rpm = DirichletProcess(alpha)

z = tzeros(Int, N)
cluster_counts = tzeros(Int, N)
fill!(cluster_counts, 0)

for i in 1:N
z[i] ~ ChineseRestaurantProcess(rpm, cluster_counts)
cluster_counts[z[i]] += 1
end

Kmax = findlast(!iszero, cluster_counts)
m = M(undef, Kmax)
for k = 1:Kmax
m[k] ~ Normal(1.0, 1.0)
end
end
model = imm(randn(100), 1.0);
sample(model, Gibbs(MH(10, :z), HMC(0.01, 4, :m)), 100);
sample(model, Gibbs(PG(10, :z), HMC(0.01, 4, :m)), 100);
end
end