In the previous post I gave a very quick introduction to the smfsb R package. As mentioned in that post, although good for teaching and learning, R isn’t a great language for serious scientific computing or computational statistics. So for the publication of the third edition of my textbook, Stochastic modelling for systems biology, I have created a library in the Scala programming language replicating the functionality provided by the R package. Here I will give a very quick introduction to the scala-smfsb library. Some familiarity with both Scala and the smfsb R package will be helpful, but is not strictly necessary. Note that the library relies on the Scala Breeze library for linear algebra and probability distributions, so some familiarity with that library can also be helpful.
Setup
To follow the along you need to have Sbt installed, and this in turn requires a recent JDK. If you are new to Scala, you may find the setup page for my Scala course to be useful, but note that on many Linux systems it can be as simple as installing the packages openjdk-8-jdk and sbt.
Once you have Sbt installed, you should be able to run it by entering sbt at your OS command line. You now need to use Sbt to create a Scala REPL with a dependency on the scala-smfsb library. There are many ways to do this, but if you are new to Scala, the simplest way is probably to start up Sbt from an empty or temporary directory (which doesn’t contain any Scala code), and then paste the following into the Sbt prompt:
set libraryDependencies += "com.github.darrenjw" %% "scala-smfsb" % "0.6" set libraryDependencies += "org.scalanlp" %% "breeze-viz" % "0.13.2" set scalaVersion := "2.12.6" set scalacOptions += "-Yrepl-class-based" console
The first time you run this it will take a little while to download and cache various library dependencies. But everything is cached, so it should be much quicker in future. When it is finished, you should have a Scala REPL ready to enter Scala code.
An introduction to scala-smfsb
It should be possible to type or copy-and-paste the commands below one-at-a-time into the Scala REPL. We need to start with a few imports.
import smfsb._
import breeze.linalg.{Vector => BVec, _}
import breeze.numerics._
import breeze.plot._
Note that I’ve renamed Breeze’s Vector type to BVec to avoid clashing with that in the Scala standard library. We are now ready to go.
Simulating models
Let’s begin by instantiating a Lotka-Volterra model, simulating a single realisation of the process, and then plotting it.
// Simulate LV with Gillespie val model = SpnModels.lv[IntState]() val step = Step.gillespie(model) val ts = Sim.ts(DenseVector(50, 100), 0.0, 20.0, 0.05, step) Sim.plotTs(ts, "Gillespie simulation of LV model with default parameters")
The library comes with a few other models. There’s a Michaelis-Menten enzyme kinetics model:
// Simulate other models with Gillespie val stepMM = Step.gillespie(SpnModels.mm[IntState]()) val tsMM = Sim.ts(DenseVector(301,120,0,0), 0.0, 100.0, 0.5, stepMM) Sim.plotTs(tsMM, "Gillespie simulation of the MM model")
and an auto-regulatory genetic network model, for example.
val stepAR = Step.gillespie(SpnModels.ar[IntState]()) val tsAR = Sim.ts(DenseVector(10, 0, 0, 0, 0), 0.0, 500.0, 0.5, stepAR) Sim.plotTs(tsAR, "Gillespie simulation of the AR model")
If you know the book and/or the R package, these models should all be familiar.
We are not restricted to exact stochastic simulation using the Gillespie algorithm. We can use an approximate Poisson time-stepping algorithm.
// Simulate LV with other algorithms val stepPts = Step.pts(model) val tsPts = Sim.ts(DenseVector(50, 100), 0.0, 20.0, 0.05, stepPts) Sim.plotTs(tsPts, "Poisson time-step simulation of the LV model")
Alternatively, we can instantiate the example models using a continuous state rather than a discrete state, and then simulate using algorithms based on continous approximations, such as Euler-Maruyama simulation of a chemical Langevin equation (CLE) approximation.
val stepCle = Step.cle(SpnModels.lv[DoubleState]()) val tsCle = Sim.ts(DenseVector(50.0, 100.0), 0.0, 20.0, 0.05, stepCle) Sim.plotTs(tsCle, "Euler-Maruyama/CLE simulation of the LV model")
If we want to ignore noise temporarily, there’s also a simple continuous deterministic Euler integrator built-in.
val stepE = Step.euler(SpnModels.lv[DoubleState]()) val tsE = Sim.ts(DenseVector(50.0, 100.0), 0.0, 20.0, 0.05, stepE) Sim.plotTs(tsE, "Continuous-deterministic Euler simulation of the LV model")
Spatial stochastic reaction-diffusion simulation
We can do 1d reaction-diffusion simulation with something like:
val N = 50; val T = 40.0 val model = SpnModels.lv[IntState]() val step = Spatial.gillespie1d(model,DenseVector(0.8, 0.8)) val x00 = DenseVector(0, 0) val x0 = DenseVector(50, 100) val xx00 = Vector.fill(N)(x00) val xx0 = xx00.updated(N/2,x0) val output = Sim.ts(xx0, 0.0, T, 0.2, step) Spatial.plotTs1d(output)
For 2d simulation, we use PMatrix, a comonadic matrix/image type defined within the library, with parallelised map and coflatMap (cobind) operations. See my post on comonads for scientific computing for further details on the concepts underpinning this, though note that it isn’t necessary to understand comonads to use the library.
val r = 20; val c = 30
val model = SpnModels.lv[DoubleState]()
val step = Spatial.cle2d(model, DenseVector(0.6, 0.6), 0.05)
val x00 = DenseVector(0.0, 0.0)
val x0 = DenseVector(50.0, 100.0)
val xx00 = PMatrix(r, c, Vector.fill(r*c)(x00))
val xx0 = xx00.updated(c/2, r/2, x0)
val output = step(xx0, 0.0, 8.0)
val f = Figure("2d LV reaction-diffusion simulation")
val p0 = f.subplot(2, 1, 0)
p0 += image(PMatrix.toBDM(output map (_.data(0))))
val p1 = f.subplot(2, 1, 1)
p1 += image(PMatrix.toBDM(output map (_.data(1))))
Bayesian parameter inference
The library also includes functions for carrying out parameter inference for stochastic dynamical systems models, using particle MCMC, ABC and ABC-SMC. See the examples directory for further details.
Next steps
Having worked through this post, the next step is to work through the tutorial. There is some overlap of content with this blog post, but the tutorial goes into more detail regarding the basics. It also finishes with suggestions for how to proceed further.
Source
This post started out as a tut document (the Scala equivalent of an RMarkdown document). The source can be found here.
of the process is being updated. That is, we target 
(for known, fixed, 
) rather than 
. This special case is known as the particle independent Metropolis-Hastings (PIMH) sampler.
using a bootstrap filter, and then accepting the proposal with probability 
, where 
is the 
is the bootstrap filter’s estimate of marginal likelihood for the new path, and 
is the estimate associated with the current path. Again using notation from the 
. Again, be sure to see the previous post for the explanation.
rather than the exact conditional distribution 
![\displaystyle \frac{\tilde{q}(\mathbf{x}_0,\ldots,\mathbf{x}_T,\mathbf{a}_0,\ldots,\mathbf{a}_{T-1})} {\displaystyle p(x_0^{b_0^k})\left[\prod_{t=0}^{T-1} \pi_t^{b_t^k}p\left(x_{t+1}^{b_{t+1}^k}|x_t^{b_t^k}\right)\right]}](https://s0.wp.com/latex.php?latex=%5Cdisplaystyle+%5Cfrac%7B%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%7D+%7B%5Cdisplaystyle+p%28x_0%5E%7Bb_0%5Ek%7D%29%5Cleft%5B%5Cprod_%7Bt%3D0%7D%5E%7BT-1%7D+%5Cpi_t%5E%7Bb_t%5Ek%7Dp%5Cleft%28x_%7Bt%2B1%7D%5E%7Bb_%7Bt%2B1%7D%5Ek%7D%7Cx_t%5E%7Bb_t%5Ek%7D%5Cright%29%5Cright%5D%7D&bg=ffffff&fg=333333&s=0&c=20201002)
.
paths, of which the conditioned-on trajectory is one.
be the path that is to be conditioned on, with ancestral lineage 
. Then, for 
, sample 
and set 
. Now suppose that at time 
we have a weighted sample from 
. First resample by sampling 
. Next sample 
. Then for all 
set 
and normalise with 
. Propagate this weighted set of particles to the next time point. At time 
select a single trajectory by sampling 
.
of the 
random variables in the augmented state space at each iteration. Since the block of variables to be updated is random, this defines an ergodic sampler for 
particles, and we have explained why the marginal distribution of the selected trajectory is the exact conditional distribution.
and the currently sampled path, 
. On the other hand, the PMMH algorithm is an MCMC algorithm which targets the full joint posterior distribution 
). 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 
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 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, 
, 
, 
. Initialise the particle filter with 
, where 
and 
is undefined). Now suppose at time 
. First resample by sampling 
, 
. Here we use 
for the discrete distribution on 
with probability mass function 
. Next sample 
. Set 
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 
. 
. 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 
Darren Wilkinson's blog 