MH Constructor

Project.toml

@@ -44,7 +44,7 @@ AdvancedMH = "0.6.8, 0.7"
 AdvancedPS = "0.4"
 AdvancedVI = "0.2"
 BangBang = "0.3"
-Bijectors = "0.12"
+Bijectors = "0.13.2"
 DataStructures = "0.18"
 Distributions = "0.23.3, 0.24, 0.25"
 DistributionsAD = "0.6"

src/Turing.jl

@@ -73,7 +73,6 @@ export  @model,                 # modelling
         Prior,                  # Sampling from the prior
 
         MH,                     # classic sampling
-        RWMH,
         Emcee,
         ESS,
         Gibbs,

src/inference/Inference.jl

@@ -99,6 +99,24 @@ Wrap a sampler so it can be used as an inference algorithm.
 """
 externalsampler(sampler::AbstractSampler) = ExternalSampler(sampler)
 
+"""
+    ESLogDensityFunction
+
+A log density function for the External sampler.
+
+"""
+const ESLogDensityFunction{M<:Model,S<:Sampler{<:ExternalSampler},V<:AbstractVarInfo} = Turing.LogDensityFunction{V,M,<:DynamicPPL.DefaultContext}
+function LogDensityProblems.logdensity(f::ESLogDensityFunction, x::NamedTuple)
+    return DynamicPPL.logjoint(f.model, DynamicPPL.unflatten(f.varinfo, x))
+end
+
+# TODO: make a nicer `set_namedtuple!` and move these functions to DynamicPPL.
+function DynamicPPL.unflatten(vi::TypedVarInfo, θ::NamedTuple) 
+    set_namedtuple!(deepcopy(vi), θ)
+    return vi
+end
+DynamicPPL.unflatten(vi::SimpleVarInfo, θ::NamedTuple) = SimpleVarInfo(θ, vi.logp, vi.transformation)
+
 # Algorithm for sampling from the prior
 struct Prior <: InferenceAlgorithm end
 

src/inference/mh.jl

@@ -188,6 +188,20 @@ function MH(space...)
     return MH{tuple(syms...), typeof(proposals)}(proposals)
 end
 
+# Some of the proposals require working in unconstrained space.
+transform_maybe(proposal::AMH.Proposal) = proposal
+function transform_maybe(proposal::AMH.RandomWalkProposal)
+    return AMH.RandomWalkProposal(Bijectors.transformed(proposal.proposal))
+end
+
+function MH(model::Model; proposal_type=AMH.StaticProposal)
+    priors = DynamicPPL.extract_priors(model)
+    props = Tuple([proposal_type(prop) for prop in values(priors)])
+    vars = Tuple(map(Symbol, collect(keys(priors))))
+    priors = map(transform_maybe, NamedTuple{vars}(props))
+    return AMH.MetropolisHastings(priors)
+end
+
 #####################
 # Utility functions #
 #####################
@@ -346,6 +360,7 @@ end
 function should_link(varinfo, sampler, proposal::AdvancedMH.RandomWalkProposal)
     return true
 end
+# FIXME: This won't be hit unless `vals` are all the exactly same concrete type of `AdvancedMH.RandomWalkProposal`!
 function should_link(
     varinfo,
     sampler,

test/contrib/inference/abstractmcmc.jl

@@ -41,7 +41,7 @@ function initialize_mh(model)
 end
 
 @testset "External samplers" begin
-    @testset "AdvancedHMC.jl" begin
+    @turing_testset "AdvancedHMC.jl" begin
         for model in DynamicPPL.TestUtils.DEMO_MODELS
             # Need some functionality to initialize the sampler.
             # TODO: Remove this once the constructors in the respective packages become "lazy".
@@ -52,12 +52,13 @@ end
                 5_000;
                 n_adapts=1_000,
                 discard_initial=1_000,
-                rtol=0.2
+                rtol=0.2,
+                sampler_name="AdvancedHMC"
             )
         end
     end
 
-    @testset "AdvancedMH.jl" begin
+    @turing_testset "AdvancedMH.jl" begin
         for model in DynamicPPL.TestUtils.DEMO_MODELS
             # Need some functionality to initialize the sampler.
             # TODO: Remove this once the constructors in the respective packages become "lazy".
@@ -68,7 +69,8 @@ end
                 10_000;
                 discard_initial=1_000,
                 thinning=10,
-                rtol=0.2
+                rtol=0.2,
+                sampler_name="AdvancedMH"
             )
         end
     end

test/inference/mh.jl

@@ -17,6 +17,12 @@
 
         s4 = Gibbs(MH(:m), MH(:s))
         c4 = sample(gdemo_default, s4, N)
+
+        s5 = externalsampler(MH(gdemo_default, proposal_type=AdvancedMH.RandomWalkProposal))
+        c5 = sample(gdemo_default, s5, N)
+
+        s6 = externalsampler(MH(gdemo_default, proposal_type=AdvancedMH.StaticProposal))
+        c6 = sample(gdemo_default, s6, N)
     end
     @numerical_testset "mh inference" begin
         Random.seed!(125)

test/modes/OptimInterface.jl

@@ -1,38 +1,3 @@
-# TODO: Remove these once the equivalent is present in `DynamicPPL.TestUtils.
-function likelihood_optima(::DynamicPPL.TestUtils.UnivariateAssumeDemoModels)
-    return (s=1/16, m=7/4)
-end
-function posterior_optima(::DynamicPPL.TestUtils.UnivariateAssumeDemoModels)
-    # TODO: Figure out exact for `s`.
-    return (s=0.907407, m=7/6)
-end
-
-function likelihood_optima(model::DynamicPPL.TestUtils.MultivariateAssumeDemoModels)
-    # Get some containers to fill.
-    vals = Random.rand(model)
-
-    # NOTE: These are "as close to zero as we can get".
-    vals.s[1] = 1e-32
-    vals.s[2] = 1e-32
-
-    vals.m[1] = 1.5
-    vals.m[2] = 2.0
-
-    return vals
-end
-function posterior_optima(model::DynamicPPL.TestUtils.MultivariateAssumeDemoModels)
-    # Get some containers to fill.
-    vals = Random.rand(model)
-
-    # TODO: Figure out exact for `s[1]`.
-    vals.s[1] = 0.890625
-    vals.s[2] = 1
-    vals.m[1] = 3/4
-    vals.m[2] = 1
-
-    return vals
-end
-
 # Used for testing how well it works with nested contexts.
 struct OverrideContext{C,T1,T2} <: DynamicPPL.AbstractContext
     context::C
@@ -57,7 +22,7 @@ function DynamicPPL.tilde_observe(context::OverrideContext, right, left, vi)
     return context.loglikelihood_weight, vi
 end
 
-@testset "OptimInterface.jl" begin
+@numerical_testset "OptimInterface.jl" begin
     @testset "MLE" begin
         Random.seed!(222)
         true_value = [0.0625, 1.75]
@@ -157,7 +122,7 @@ end
     # FIXME: Some models doesn't work for Tracker and ReverseDiff.
     if Turing.Essential.ADBACKEND[] === :forwarddiff
         @testset "MAP for $(model.f)" for model in DynamicPPL.TestUtils.DEMO_MODELS
-            result_true = posterior_optima(model)
+            result_true = DynamicPPL.TestUtils.posterior_optima(model)
 
             @testset "$(nameof(typeof(optimizer)))" for optimizer in [LBFGS(), NelderMead()]
                 result = optimize(model, MAP(), optimizer)
@@ -188,7 +153,7 @@ end
             DynamicPPL.TestUtils.demo_dot_assume_matrix_dot_observe_matrix,
         ]
         @testset "MLE for $(model.f)" for model in DynamicPPL.TestUtils.DEMO_MODELS
-            result_true = likelihood_optima(model)
+            result_true = DynamicPPL.TestUtils.likelihood_optima(model)
 
             # `NelderMead` seems to struggle with convergence here, so we exclude it.
             @testset "$(nameof(typeof(optimizer)))" for optimizer in [LBFGS(),]