diff --git a/src/mh-core.jl b/src/mh-core.jl index cfeb24e..abdf03b 100644 --- a/src/mh-core.jl +++ b/src/mh-core.jl @@ -206,8 +206,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α @@ -216,3 +215,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 2a4a05d..7b90435 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -44,6 +44,10 @@ using ForwardDiff 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=["μ", "σ"])