Introduction
In the previous post I gave a brief introduction to Rainier, a new HMC-based probabilistic programming library/DSL for Scala. In that post I assumed that people were using the latest source version of the library. Since then, version 0.1.1 of the library has been released, so in this post I will demonstrate use of the released version of the software (using the binaries published to Sonatype), and will walk through a slightly more interesting example – a dynamic linear state space model with unknown static parameters. This is similar to, but slightly different from, the DLM example in the Rainier library. So to follow along with this post, all that is required is SBT.
An interactive session
First run SBT from an empty directory, and paste the following at the SBT prompt:
set libraryDependencies += "com.stripe" %% "rainier-plot" % "0.1.1" set scalaVersion := "2.12.4" console
This should give a Scala REPL with appropriate dependencies (rainier-plot has all of the relevant transitive dependencies). We’ll begin with some imports, and then simulating some synthetic data from a dynamic linear state space model with an AR(1) latent state and Gaussian noise on the observations.
import com.stripe.rainier.compute._ import com.stripe.rainier.core._ import com.stripe.rainier.sampler._ implicit val rng = ScalaRNG(1) val n = 60 // number of observations/time points val mu = 3.0 // AR(1) mean val a = 0.95 // auto-regressive parameter val sig = 0.2 // AR(1) SD val sigD = 3.0 // observational SD val state = Stream. iterate(0.0)(x => mu + (x - mu) * a + sig * rng.standardNormal). take(n).toVector val obs = state.map(_ + sigD * rng.standardNormal)
Now we have some synthetic data, let’s think about building a probabilistic program for this model. Start with a prior.
case class Static(mu: Real, a: Real, sig: Real, sigD: Real)
val prior = for {
mu <- Normal(0, 10).param
a <- Normal(1, 0.1).param
sig <- Gamma(2,1).param
sigD <- Gamma(2,2).param
sp <- Normal(0, 50).param
} yield (Static(mu, a, sig, sigD), List(sp))
Note the use of a case class for wrapping the static parameters. Next, let’s define a function to add a state and associated observation to an existing model.
def addTimePoint(current: RandomVariable[(Static, List[Real])],
datum: Double) = for {
tup <- current
static = tup._1
states = tup._2
os = states.head
ns <- Normal(((Real.one - static.a) * static.mu) + (static.a * os),
static.sig).param
_ <- Normal(ns, static.sigD).fit(datum)
} yield (static, ns :: states)
Given this, we can generate the probabilistic program for our model as a fold over the data initialised with the prior.
val fullModel = obs.foldLeft(prior)(addTimePoint(_, _))
If we don’t want to keep samples for all of the variables, we can focus on the parameters of interest, wrapping the results in a Map for convenient sampling and plotting.
val model = for {
tup <- fullModel
static = tup._1
states = tup._2
} yield
Map("mu" -> static.mu,
"a" -> static.a,
"sig" -> static.sig,
"sigD" -> static.sigD,
"SP" -> states.reverse.head)
We can sample with
val out = model.sample(HMC(3), 100000, 10000 * 500, 500)
(this will take several minutes) and plot some diagnostics with
import com.cibo.evilplot.geometry.Extent
import com.stripe.rainier.plot.EvilTracePlot._
val truth = Map("mu" -> mu, "a" -> a, "sigD" -> sigD,
"sig" -> sig, "SP" -> state(0))
render(traces(out, truth), "traceplots.png",
Extent(1200, 1400))
render(pairs(out, truth), "pairs.png")
This generates the following diagnostic plots:
Everything looks good.
Summary
Rainier is a monadic embedded DSL for probabilistic programming in Scala. We can use standard functional combinators and for-expressions for building models to sample, and then run an efficient HMC algorithm on the resulting probability monad in order to obtain samples from the posterior distribution of the model.
See the Rainier repo for further details.
random quantities to construct a Monte Carlo estimate of the integral
.
separate computations, and use a map to compute the required Monte Carlo sums.
, latent variables 
and data 
, of the form
. On the other hand, the PMMH algorithm is an MCMC algorithm which targets the full joint posterior distribution 
. Now, the PMMH scheme does reduce to the pseudo-marginal scheme if samples of 
). Below I will describe the algorithm and explain why it works, but first it is necessary to understand the relationship between marginal, joint and “likelihood-free” MCMC updating schemes for such latent variable models.
is proposed from 
and accepted with probability 
, where
required in the acceptance ratio is often difficult to compute.
even when we can’t evaluate it, or marginalise 
is simulated from the model 
. The pair 
is then jointly accepted with ratio
. The pair 
, as the algorithms above. First, a proposed new 
is generated by running a bootstrap particle filter (as described in the previous post, and below) using the proposed new model parameters, 
is accepted using the Metropolis-Hastings ratio
is the particle filter’s (unbiased) estimate of marginal likelihood, described in the previous post, and below. Note that this approach tends to the perfect joint/marginal updating scheme as the number of particles used in the filter tends to infinity. Note also that for a single particle, the particle filter just blindly forward simulates from 
and that the filter’s estimate of marginal likelihood is just the observed data likelihood 
leading precisely to the simple likelihood-free scheme. To understand for an arbitrary finite number of particles, 
, one needs to think carefully about the structure of the particle filter.
, 
, 
. Initialise the particle filter with 
, where 
and 
(note that 
is undefined). Now suppose at time 
we have a sample from 
: 
. First resample by sampling 
, 
. Here we use 
for the discrete distribution on 
with probability mass function 
. Next sample 
. Set 
and 
. Finally, propagate 
to the next step… We define the filter’s estimate of likelihood as 
and 
. See ![\displaystyle \tilde{q}(\mathbf{x}_0,\ldots,\mathbf{x}_T,\mathbf{a}_0,\ldots,\mathbf{a}_{T-1}) = \left[\prod_{k=1}^M p(x_0^k)\right] \left[\prod_{t=0}^{T-1} \prod_{k=1}^M \pi_t^{a_t^k} p(x_{t+1}^k|x_t^{a_t^k}) \right]](https://s0.wp.com/latex.php?latex=%5Cdisplaystyle++%5Ctilde%7Bq%7D%28%5Cmathbf%7Bx%7D_0%2C%5Cldots%2C%5Cmathbf%7Bx%7D_T%2C%5Cmathbf%7Ba%7D_0%2C%5Cldots%2C%5Cmathbf%7Ba%7D_%7BT-1%7D%29++%3D+%5Cleft%5B%5Cprod_%7Bk%3D1%7D%5EM+p%28x_0%5Ek%29%5Cright%5D+%5Cleft%5B%5Cprod_%7Bt%3D0%7D%5E%7BT-1%7D++++%5Cprod_%7Bk%3D1%7D%5EM+%5Cpi_t%5E%7Ba_t%5Ek%7D+p%28x_%7Bt%2B1%7D%5Ek%7Cx_t%5E%7Ba_t%5Ek%7D%29+%5Cright%5D++&bg=ffffff&fg=333333&s=0&c=20201002)
from 
giving the joint density 
, and 
. If we now think about the structure of the PMMH algorithm, our proposal on the space of all random variables in the problem is in fact
when we consider just the parameters and the selected trajectory. But if we consider the terms in the joint distribution of the proposal corresponding to the trajectory selected by ![\displaystyle p_\theta(x_0^{b_0^{k'}})\left[\prod_{t=0}^{T-1} \pi_t^{b_t^{k'}} p_\theta(x_{t+1}^{b_{t+1}^{k'}}|x_t^{b_t^{k'}})\right]\pi_T^{k'} = p_\theta(x_{0:T}^{k'})\prod_{t=0}^T \pi_t^{b_t^{k'}}](https://s0.wp.com/latex.php?latex=%5Cdisplaystyle++p_%5Ctheta%28x_0%5E%7Bb_0%5E%7Bk%27%7D%7D%29%5Cleft%5B%5Cprod_%7Bt%3D0%7D%5E%7BT-1%7D+%5Cpi_t%5E%7Bb_t%5E%7Bk%27%7D%7D++++p_%5Ctheta%28x_%7Bt%2B1%7D%5E%7Bb_%7Bt%2B1%7D%5E%7Bk%27%7D%7D%7Cx_t%5E%7Bb_t%5E%7Bk%27%7D%7D%29%5Cright%5D%5Cpi_T%5E%7Bk%27%7D++%3D++p_%5Ctheta%28x_%7B0%3AT%7D%5E%7Bk%27%7D%29%5Cprod_%7Bt%3D0%7D%5ET+%5Cpi_t%5E%7Bb_t%5E%7Bk%27%7D%7D++&bg=ffffff&fg=333333&s=0&c=20201002)
in terms of the unnormalised weights, simplifies to
in the acceptance ratio, we knock out the normalising constant, allowing all of the other terms in the proposal to be marginalised away. In other words, the target of the chain can be written as proportional to
together with a prior on 
, giving a joint model 
from some fairly arbitrary proposal distribution and then accepting the value with probability 
will still lead to the exact posterior if the estimate is unbiased in the sense that ![E[\hat\pi(y|\theta)]=\pi(y|\theta)](https://s0.wp.com/latex.php?latex=E%5B%5Chat%5Cpi%28y%7C%5Ctheta%29%5D%3D%5Cpi%28y%7C%5Ctheta%29&bg=ffffff&fg=333333&s=0&c=20201002)
. Consequently, sources of (cheap) unbiased Monte Carlo estimates of (marginal) likelihood are of potential interest in the development of exact MCMC algorithms.
and simulate realisations from 
. That is, simulate values 
from 
, and then put
. That is, 
having the same support as 
. The weights, 
, are known as importance weights.
using values sampled from an auxiliary distribution 
, where we now supress any dependence of the distributions on model parameters, 
from 
. Then compute normalised weights 
. Generate a new sample of size 
case, the final sample will be exactly drawn from the auxiliary and not the target. The procedure is asymptotic, in that it improves as the sample size increases, tending to the exact target as 
. The proportion of the auxiliary samples falling in a small interval 
will be 
, corresponding to roughly 
particles. The weight for each of those particles will be 
, and since the expected weight of a random particle is 1, the sum of all weights will be (roughly) ![\tilde{w}(x)=\pi(x)/[N\pi'(x)]](https://s0.wp.com/latex.php?latex=%5Ctilde%7Bw%7D%28x%29%3D%5Cpi%28x%29%2F%5BN%5Cpi%27%28x%29%5D&bg=ffffff&fg=333333&s=0&c=20201002)
. The combined weight of all particles in 
. Clearly then, when we resample 
particles from this interval. This corresponds to a proportion 
, corresponding to a density of 
governed by a transition kernel 
is partially observed via some measurement model 
leading to data 
. The idea is to make inference for the hidden states 
. The method is a very simple application of the importance resampling technique. At each time, 
.
with uniform normalised weights 
. Then suppose that we have a weighted sample 
from 
(giving an approximate random sample from 
. Next propagate each particle forward according to the Markov process model by sampling 
(giving an approximate random sample from 
). Then for each of the new particles, compute a weight 
. 
is a consistent estimator of 
. It is much less clear, but nevertheless true that this estimator is also unbiased. The standard reference for this fact is 
. This turns out to be a simple special case of the particle marginal Metropolis-Hastings (PMMH) algorithm described in 
using a proposal 
. As explained previously, the required Metropolis-Hastings acceptance ratio will take the form 
, can be computed. We can obviously just substitute this estimate into the acceptance ratio to get 
, where the expectation is with respect to the Monte-Carlo error for a given fixed value of 
, representing the noise in the Monte-Carlo estimate of 
(note that in an abuse of notation, the function 
is unrelated to 
, where 
is a constant independent of 
, we have the (typical) special case of 
, but we will consider relaxing this constraint later. 
is being proposed in addition to a new value for 
, it is clear that the proposal generates 
from the density 
, and we can re-write our “approximate” acceptance ratio as 
, but the marginal for 
Darren Wilkinson's blog 