From fe1091d0b4391a5af273a3dc25ea533e54fddc5c Mon Sep 17 00:00:00 2001 From: Arthur Lui Date: Tue, 8 Dec 2020 13:40:56 -0700 Subject: [PATCH] Attempt to fix #41. --- src/mh-core.jl | 66 ++++++++++++++++++++++++++++++++++++++++++++++-- test/runtests.jl | 4 +++ 2 files changed, 68 insertions(+), 2 deletions(-) diff --git a/src/mh-core.jl b/src/mh-core.jl index 33832a2..fb74f09 100644 --- a/src/mh-core.jl +++ b/src/mh-core.jl @@ -203,8 +203,7 @@ function AbstractMCMC.step( params = propose(rng, spl, model, params_prev) # Calculate the log acceptance probability. - logα = logdensity(model, params) - logdensity(model, params_prev) + - q(spl, params_prev, params) - q(spl, params, params_prev) + logα = compute_log_acceptance_prob(model, params, spl, params_prev) # Decide whether to return the previous params or the new one. if -Random.randexp(rng) < logα @@ -213,3 +212,66 @@ function AbstractMCMC.step( return params_prev, params_prev end end + +""" +Type alias for symmetric distributions. This is not an exhaustive list. (e.g. +`Uniform` can be symmetric.) +""" +const SymmetricDistribution = Union{Normal,MvNormal,MvTDist,Cauchy,TDist} + +""" +Type alias for symmetric random walk proposals +""" +const SymmetricRandomWalkProposal = let + Union{RandomWalkProposal{<:SymmetricDistribution}, + RandomWalkProposal{<:AbstractVector{<:SymmetricDistribution}}} +end + +""" +Type alias for symmetric proposals. + +NOTE to developers: Currently, this is `SymmetricRandomWalkProposal`, but new +proposals (e.g. an adaptive random walk proposal) can be added to the `Union`. +""" +const SymmetricProposal = Union{SymmetricRandomWalkProposal} + +""" +Computes metropolis acceptance ratio (without Hastings). +""" +function compute_log_metropolis_ratio( + model::DensityModel, + params::T, + spl::MetropolisHastings, + params_prev::Transition +) where T + return logdensity(model, params) - logdensity(model, params_prev) +end + +""" +Computes log acceptance ratio for symmetric proposals. This is simply +the log Metropolis ratio. +""" +function compute_log_acceptance_prob( + model::DensityModel, + params::T, + spl::MetropolisHastings{<:SymmetricProposal}, + params_prev::Transition) where T + + # Don't compute Hastings ratio for symmetric proposal distributions. + return compute_log_metropolis_ratio(model, params, spl, params_prev) +end + +""" +Computes log acceptance ratio for symmetric proposals. This is the +full log-Metropolis-Hastings acceptance probability. +""" +function compute_log_acceptance_prob( + model::DensityModel, + params::T, + spl::MetropolisHastings, + params_prev::Transition) where T + + metropolis_ratio = compute_log_metropolis_ratio(model, params, spl, params_prev) + hastings_ratio = q(spl, params_prev, params) - q(spl, params, params_prev) + return metropolis_ratio + hastings_ratio +end diff --git a/test/runtests.jl b/test/runtests.jl index 2e81db0..639513a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -42,6 +42,10 @@ using Test spl1 = RWMH([Normal(0,1), Normal(0, 1)]) spl2 = RWMH(MvNormal([0.0, 0.0], 1)) + # Ensure that these are `SymmetricProposal`s. + @test spl1 isa AdvancedMH.MetropolisHastings{<:AdvancedMH.SymmetricProposal} + @test spl2 isa AdvancedMH.MetropolisHastings{<:AdvancedMH.SymmetricProposal} + # Sample from the posterior. chain1 = sample(model, spl1, 100000; chain_type=StructArray, param_names=["μ", "σ"]) chain2 = sample(model, spl2, 100000; chain_type=StructArray, param_names=["μ", "σ"])