Merge pull request #23 from devmotion/suggestions

Some suggestions

Author: cpfiffer

Project.toml

@@ -1,6 +1,6 @@
 name = "AdvancedMH"
 uuid = "5b7e9947-ddc0-4b3f-9b55-0d8042f74170"
-version = "0.4.1"
+version = "0.5.0"
 
 [deps]
 AbstractMCMC = "80f14c24-f653-4e6a-9b94-39d6b0f70001"

README.md

@@ -75,34 +75,35 @@ supported methods are `Array{Proposal}`, `Proposal`, and `NamedTuple{Proposal}`.
 ```julia
 # Provide a univariate proposal.
 m1 = DensityModel(x -> logpdf(Normal(x,1), 1.0))
-p1 = Proposal(Static(), Normal(0,1))
+p1 = StaticProposal(Normal(0,1))
 c1 = sample(m1, MetropolisHastings(p1), 100; chain_type=Vector{NamedTuple})
 
 # Draw from a vector of distributions.
 m2 = DensityModel(x -> logpdf(Normal(x[1], x[2]), 1.0))
-p2 = Proposal(Static(), [Normal(0,1), InverseGamma(2,3)])
+p2 = StaticProposal([Normal(0,1), InverseGamma(2,3)])
 c2 = sample(m2, MetropolisHastings(p2), 100; chain_type=Vector{NamedTuple})
 
 # Draw from a `NamedTuple` of distributions.
 m3 = DensityModel(x -> logpdf(Normal(x.a, x.b), 1.0))
-p3 = (a=Proposal(Static(), Normal(0,1)), b=Proposal(Static(), InverseGamma(2,3)))
+p3 = (a=StaticProposal(Normal(0,1)), b=StaticProposal(InverseGamma(2,3)))
 c3 = sample(m3, MetropolisHastings(p3), 100; chain_type=Vector{NamedTuple})
 
 # Draw from a functional proposal.
 m4 = DensityModel(x -> logpdf(Normal(x,1), 1.0))
-p4 = Proposal(Static(), (x=1.0) -> Normal(x, 1))
+p4 = StaticProposal((x=1.0) -> Normal(x, 1))
 c4 = sample(m4, MetropolisHastings(p4), 100; chain_type=Vector{NamedTuple})
 ```
 
 ## Static vs. Random Walk
 
 Currently there are only two methods of inference available. Static MH simply draws from the prior, with no
 conditioning on the previous sample. Random walk will add the proposal to the previously observed value.
-If you are constructing a `Proposal` by hand, you can determine whether the proposal is `Static` or `RandomWalk` using
+If you are constructing a `Proposal` by hand, you can determine whether the proposal is a
+`StaticProposal` or a `RandomWalkProposal` using
 
 ```julia
-static_prop = Proposal(Static(), Normal(0,1))
-rw_prop = Proposal(RandomWalk(), Normal(0,1))
+static_prop = StaticProposal(Normal(0,1))
+rw_prop = RandomWalkProposal(Normal(0,1))
 ```
 
 Different methods are easily composeable. One parameter can be static and another can be a random walk,

src/AdvancedMH.jl

@@ -5,21 +5,17 @@ import AbstractMCMC
 using Distributions
 using Requires
 
-using Random
+import Random
 
 # Exports
-export MetropolisHastings, DensityModel, RWMH, StaticMH, Proposal, Static, RandomWalk
+export MetropolisHastings, DensityModel, RWMH, StaticMH, StaticProposal, RandomWalkProposal
 
 # Reexports
 using AbstractMCMC: sample, psample
 export sample, psample
 
 # Abstract type for MH-style samplers.
 abstract type Metropolis <: AbstractMCMC.AbstractSampler end
-abstract type ProposalStyle end
-
-struct RandomWalk <: ProposalStyle end
-struct Static <: ProposalStyle end
 
 # Define a model type. Stores the log density function and the data to 
 # evaluate the log density on.
@@ -55,7 +51,7 @@ logdensity(model::DensityModel, t::Transition) = t.lp
 
 # A basic chains constructor that works with the Transition struct we defined.
 function AbstractMCMC.bundle_samples(
-    rng::AbstractRNG, 
+    rng::Random.AbstractRNG, 
     model::DensityModel, 
     s::Metropolis, 
     N::Integer,
@@ -68,7 +64,7 @@ function AbstractMCMC.bundle_samples(
 end
 
 function AbstractMCMC.bundle_samples(
-    rng::AbstractRNG, 
+    rng::Random.AbstractRNG, 
     model::DensityModel, 
     s::Metropolis, 
     N::Integer, 

src/mcmcchains-connect.jl

@@ -2,7 +2,7 @@ import .MCMCChains: Chains
 
 # A basic chains constructor that works with the Transition struct we defined.
 function AbstractMCMC.bundle_samples(
-    rng::AbstractRNG, 
+    rng::Random.AbstractRNG, 
     model::DensityModel, 
     s::Metropolis, 
     N::Integer, 

src/mh-core.jl

@@ -12,16 +12,16 @@ x = (a = 1.0, b=3.8)
 The proposal would be
 
 ```julia
-proposal = (a=Proposal(Static(), Normal(0,1)), b=Proposal(Static(), Normal(0,1)))
+proposal = (a=StaticProposal(Normal(0,1)), b=StaticProposal(Normal(0,1)))
 ````
 
-Other allowed proposal styles are
+Other allowed proposals are
 
 ```
-p1 = Proposal(Static(), Normal(0,1))
-p2 = Proposal(Static(), [Normal(0,1), InverseGamma(2,3)])
-p3 = Proposal(Static(), (a=Normal(0,1), b=InverseGamma(2,3)))
-p4 = Proposal(Static(), (x=1.0) -> Normal(x, 1))
+p1 = StaticProposal(Normal(0,1))
+p2 = StaticProposal([Normal(0,1), InverseGamma(2,3)])
+p3 = StaticProposal(a=Normal(0,1), b=InverseGamma(2,3))
+p4 = StaticProposal((x=1.0) -> Normal(x, 1))
 ```
 
 The sampler is constructed using
@@ -41,98 +41,92 @@ used if `chain_type=Chains`.
 types are `chain_type=Chains` if `MCMCChains` is imported, or 
 `chain_type=StructArray` if `StructArrays` is imported.
 """
-mutable struct MetropolisHastings{D} <: Metropolis
-    proposal :: D
+struct MetropolisHastings{D} <: Metropolis
+    proposal::D
 end
 
-StaticMH(d) = MetropolisHastings(Proposal(Static(), d))
-RWMH(d) = MetropolisHastings(Proposal(RandomWalk(), d))
+StaticMH(d) = MetropolisHastings(StaticProposal(d))
+RWMH(d) = MetropolisHastings(RandomWalkProposal(d))
+
+# default function without RNG
+propose(spl::MetropolisHastings, args...) = propose(Random.GLOBAL_RNG, spl, args...)
 
 # Propose from a vector of proposals
 function propose(
+    rng::Random.AbstractRNG,
     spl::MetropolisHastings{<:AbstractArray},
     model::DensityModel
 )
-    proposal = map(i -> propose(spl.proposal[i], model), 1:length(spl.proposal))
+    proposal = map(p -> propose(rng, p, model), spl.proposal)
     return Transition(model, proposal)
 end
 
 function propose(
+    rng::Random.AbstractRNG,
     spl::MetropolisHastings{<:AbstractArray},
     model::DensityModel,
     params_prev::Transition
 )
-    proposal = map(i -> propose(spl.proposal[i], model, params_prev.params[i]), 1:length(spl.proposal))
+    proposal = map(spl.proposal, params_prev.params) do p, params
+        propose(rng, p, model, params)
+    end
     return Transition(model, proposal)
 end
 
 # Make a proposal from one Proposal struct.
 function propose(
+    rng::Random.AbstractRNG,
     spl::MetropolisHastings{<:Proposal},
     model::DensityModel
 )
-    proposal = propose(spl.proposal, model)
+    proposal = propose(rng, spl.proposal, model)
     return Transition(model, proposal)
 end
 
 function propose(
+    rng::Random.AbstractRNG,
     spl::MetropolisHastings{<:Proposal},
     model::DensityModel,
     params_prev::Transition
 )
-    proposal = propose(spl.proposal, model, params_prev.params)
+    proposal = propose(rng, spl.proposal, model, params_prev.params)
     return Transition(model, proposal)
 end
 
 # Make a proposal from a NamedTuple of Proposal.
 function propose(
+    rng::Random.AbstractRNG,
     spl::MetropolisHastings{<:NamedTuple},
     model::DensityModel
 )
-    proposal = _propose(spl.proposal, model)
-    return Transition(model, proposal)
-end
-
-@generated function _propose(
-    proposals::NamedTuple{names},
-    model::DensityModel
-) where {names}
-    expr = Expr(:tuple)
-    map(names) do f
-        push!(expr.args, Expr(:(=), f, :(propose(proposals.$f, model)) ))
+    proposal = map(spl.proposal) do p
+        propose(rng, p, model)
     end
-    return expr
+    return Transition(model, proposal)
 end
 
 function propose(
+    rng::Random.AbstractRNG,
     spl::MetropolisHastings{<:NamedTuple},
     model::DensityModel,
     params_prev::Transition
 )
-    proposal = _propose(spl.proposal, model, params_prev.params)
-    return Transition(model, proposal)
-end
-
-@generated function _propose(
-    proposals::NamedTuple{names},
-    model::DensityModel, 
-    params_prev::NamedTuple
-) where {names}
-    expr = Expr(:tuple)
-    map(names) do f
-        push!(expr.args, Expr(:(=), f, :(propose(proposals.$f, model, params_prev.$f)) ))
+    proposal = map(spl.proposal, params_prev.params) do p, params
+        propose(rng, p, model, params)
     end
-    return expr
+    return Transition(model, proposal)
 end
 
-
 # Evaluate the likelihood of t conditional on t_cond.
 function q(
     spl::MetropolisHastings{<:AbstractArray},
     t::Transition,
     t_cond::Transition
 )
-    return sum(map(i -> q(spl.proposal[i], t.params[i], t_cond.params[i]), 1:length(spl.proposal)))
+    # mapreduce with multiple iterators requires Julia 1.2 or later
+    return mapreduce(+, 1:length(spl.proposal)) do i
+        q(spl.proposal[i], t.params[i], t_cond.params[i])
+    end
 end
 
 function q(
@@ -148,21 +142,17 @@ function q(
     t::Transition,
     t_cond::Transition
 )
-    ks = keys(t.params)
-    total = 0.0
-
-    for k in ks
-        total += q(spl.proposal[k], t.params[k], t_cond.params[k])
+    # mapreduce with multiple iterators requires Julia 1.2 or later
+    return mapreduce(+, keys(t.params)) do k
+        q(spl.proposal[k], t.params[k], t_cond.params[k])
     end
-
-    return total
 end
 
 # Define the first step! function, which is called at the 
 # beginning of sampling. Return the initial parameter used
 # to define the sampler.
 function AbstractMCMC.step!(
-    rng::AbstractRNG,
+    rng::Random.AbstractRNG,
     model::DensityModel,
     spl::MetropolisHastings,
     N::Integer,
@@ -171,7 +161,7 @@ function AbstractMCMC.step!(
     kwargs...
 )
     if init_params === nothing
-        return propose(spl, model)
+        return propose(rng, spl, model)
     else
         return Transition(model, init_params)
     end
@@ -181,22 +171,22 @@ end
 # either a new proposal (if accepted) or the previous proposal 
 # (if not accepted).
 function AbstractMCMC.step!(
-    rng::AbstractRNG,
+    rng::Random.AbstractRNG,
     model::DensityModel,
     spl::MetropolisHastings,
     ::Integer,
     params_prev::Transition;
     kwargs...
 )
     # Generate a new proposal.
-    params = propose(spl, model, params_prev)
+    params = propose(rng, spl, model, params_prev)
 
     # Calculate the log acceptance probability.
-    α = logdensity(model, params) - logdensity(model, params_prev) + 
+    logα = logdensity(model, params) - logdensity(model, params_prev) + 
         q(spl, params_prev, params) - q(spl, params, params_prev)
 
     # Decide whether to return the previous params or the new one.
-    if log(rand(rng)) < min(0.0, α)
+    if -Random.randexp(rng) < logα
         return params
     else
         return params_prev