Reference: R.J. Boys, D.J. Wilkinson and T.B.L. Kirkwood. Bayesian inference for a discretely observed stochastic kinetic model. Statistics and Computing (2008) 18:125-135.

Statistical model

Let \(\delta_t=1/m\) where \(m\) is an integer, and \(T\) a multiple of \(m\). For \(t=1,\ldots,T\) \[ x_t|x_{t-1}\sim \mathcal N(x_{t-1}+\alpha(x_{t-1},c)\delta_t,\beta(x_{t-1},c)\delta_t) \]

where \[ \alpha(x,c) = \left( \begin{array}{c} c_1x_1-c_2x_1x_2 \\ c_2x_1x_2-c_3x_2 \\ \end{array} \right) \] and \[ \beta(x,c) = \left( \begin{array}{cc} c_1x_1+c_2x_1x_2 & -c_2x_1x_2 \\ -c_2x_1x_2 & c_2x_1x_2 + c_3x_2 \\ \end{array} \right) \]

For \(t=m,2m,3m,\ldots,T\), \[ y_t|x_t\sim \mathcal N(x_{1t},\sigma^2) \]

and for \(i=1,\ldots,3\)

\[ \log(c_i)\sim Unif(-7,2) \]

\(x_{t1}\) and \(x_{t2}\) respectively correspond to the number of preys and predators and \(y_t\) is the approximated number of preys. The model is the approximation of the Lotka-Volterra model.

Statistical model in BUGS language

model_file = 'stoch_kinetic_cle.bug' # BUGS model filename
cat(readLines(model_file), sep = "\n")
# Stochastic kinetic predator-prey model
# with chemical Langevin equations
#
# Reference: A. Golightly and D. J. Wilkinson. Bayesian parameter inference
# for stochastic biochemical network models using particle Markov chain
# Monte Carlo. Interface Focus, vol.1, pp. 807-820, 2011.

var x_true[2,t_max/dt], x_true_temp[2,t_max/dt],
    x[2,t_max/dt], x_temp[2,t_max/dt], y[t_max/dt], 
    beta[2,2,t_max/dt], beta_true[2,2,t_max/dt], 
    logc[3], c[3], c_true[3]

data
{
  x_true[,1] ~ dmnormvar(x_init_mean, x_init_var)
  for (t in 2:t_max/dt)
  {
    alpha_true[1,t] <- c_true[1]*x_true[1,t-1] - c_true[2]*x_true[1,t-1]*x_true[2,t-1]
    alpha_true[2,t] <- c_true[2]*x_true[1,t-1]*x_true[2,t-1] - c_true[3]*x_true[2,t-1]
    beta_true[1,1,t] <- c_true[1]*x_true[1,t-1] + c_true[2]*x_true[1,t-1]*x_true[2,t-1]
    beta_true[1,2,t] <- -c_true[2]*x_true[1,t-1]*x_true[2,t-1]
    beta_true[2,1,t] <- beta_true[1,2,t]
    beta_true[2,2,t] <- c_true[2]*x_true[1,t-1]*x_true[2,t-1] + c_true[3]*x_true[2,t-1]
    x_true_temp[,t] ~ dmnormvar(x_true[,t-1]+alpha_true[,t]*dt, (beta_true[,,t])*dt)
    # To avoid extinction
    x_true[1,t] <- max(x_true_temp[1,t],1)
    x_true[2,t] <- max(x_true_temp[2,t],1)
  }
  for (t in 1:t_max)
  {
    y[t/dt] ~ dnorm(x_true[1,t/dt], prec_y)
  }
}

model
{
  logc[1] ~ dunif(-7,2)
  logc[2] ~ dunif(-7,2)
  logc[3] ~ dunif(-7,2)
  c[1] <- exp(logc[1])
  c[2] <- exp(logc[2])
  c[3] <- exp(logc[3])
  x[,1] ~ dmnormvar(x_init_mean,  x_init_var)
  for (t in 2:t_max/dt)
  {
    alpha[1,t] <- c[1]*x[1,t-1] - c[2]*x[1,t-1]*x[2,t-1]
    alpha[2,t] <- c[2]*x[1,t-1]*x[2,t-1] - c[3]*x[2,t-1]
    beta[1,1,t] <- c[1]*x[1,t-1] + c[2]*x[1,t-1]*x[2,t-1]
    beta[1,2,t] <- -c[2]*x[1,t-1]*x[2,t-1]
    beta[2,1,t] <- beta[1,2,t]
    beta[2,2,t] <- c[2]*x[1,t-1]*x[2,t-1] + c[3]*x[2,t-1]
    x_temp[,t] ~ dmnormvar(x[,t-1]+alpha[,t]*dt, beta[,,t]*dt)
    # To avoid extinction
    x[1,t] <- max(x_temp[1,t],1)
    x[2,t] <- max(x_temp[2,t],1)
  }
  for (t in 1:t_max)
  {
    y[t/dt] ~ dnorm(x[1,t/dt], prec_y)
  }
}
  1. Load rbiips package
require(rbiips)

General settings

par(bty='l')
light_blue = rgb(.7, .7, 1)
light_red = rgb(1, .7, .7)
dark_blue = rgb(0, 0, .5)
dark_red = rgb(.5, 0, 0)

Set the random numbers generator seed for reproducibility

set.seed(0)

Load model and data

Model parameters

t_max = 20
dt = .2
x_init_mean = c(100, 100)
x_init_var = diag(10, 2)
c_true = c(.5, .0025, .3)
prec_y = 1/10
data = list(t_max=t_max, dt=dt, c_true=c_true,
            x_init_mean=x_init_mean, x_init_var=x_init_var, prec_y=prec_y)

Compile BUGS model and sample data

sample_data = TRUE # Boolean
model = biips_model(model_file, data, sample_data=sample_data) # Create Biips model and sample data
* Parsing model in: stoch_kinetic_cle.bug
* Compiling data graph
  Declaring variables
  Resolving undeclared variables
  Allocating nodes
  Graph size: 1914
  Sampling data
  Reading data back into data table
* Compiling model graph
  Declaring variables
  Resolving undeclared variables
  Allocating nodes
  Graph size: 2713
data = model$data()

Plot data

t_vec <- seq(dt, t_max, dt)
plot(t_vec, data$x_true[1,],
     type='l', col='blue', lwd=2,
     xlab='Time',
     ylab='Number of individuals',
     ylim=c(50,350))
lines(t_vec, data$x_true[2,], col='red', lwd=2)
points(t_vec, data$y, pch=8, col='green', lwd=2)
legend('topright', leg=c('Prey', 'Predator', 'Measurements'),
       col=c('blue', 'red', 'green'), lwd=2, lty=c(1,1,NA), pch=c(NA,NA,8),
       bty='n')
Data

Data

Biips Sensitivity analysis with Sequential Monte Carlo

Parameters of the algorithm

n_part = 100 # Number of particles
n_grid = 20
param_values = list('logc[1]' = seq(-7, 1, len=n_grid),
                    'logc[2]' = rep(log(c_true[2]), n_grid),
                    'logc[3]' = rep(log(c_true[3]), n_grid)) # Range of values

Run sensitivity analysis with SMC

out_sens = biips_smc_sensitivity(model, param_values, n_part)
* Analyzing sensitivity with 100 particles
  |--------------------------------------------------| 100%
  |**************************************************| 20 iterations in 2.99 s

Plot penalized log-marginal likelihood

plot(param_values[[1]], out_sens$log_marg_like_pen, col='blue', pch=20,
     xlab=expression(log(c[1])),
     ylab='Penalized log-marginal likelihood',
     ylim=c(-15000, 0))
Sensitivity: Penalized log-marginal likelihood

Sensitivity: Penalized log-marginal likelihood

Biips Particle Marginal Metropolis-Hastings

We now use Biips to run a Particle Marginal Metropolis-Hastings in order to obtain posterior MCMC samples of the parameters and variables \(x\).

Parameters of the PMMH

n_burn = 2000 # nb of burn-in/adaptation iterations
n_iter = 20000 # nb of iterations after burn-in
thin = 10 # thinning of MCMC outputs
n_part = 100 # nb of particles for the SMC
param_names = c('logc[1]', 'logc[2]', 'logc[3]') # names of the variables updated with MCMC (others are updated with SMC)
latent_names = c('x') # names of the variables updated with SMC and that need to be monitored

Init PMMH

obj_pmmh = biips_pmmh_init(model, param_names, inits=list(-1, -5, -1),
                           latent_names=latent_names) # creates a pmmh object
* Initializing PMMH

Run PMMH

biips_pmmh_update(obj_pmmh, n_burn, n_part) # adaptation and burn-in iterations
* Adapting PMMH with 100 particles
  |--------------------------------------------------| 100%
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 2000 iterations in 289.05 s
out_pmmh = biips_pmmh_samples(obj_pmmh, n_iter, n_part, thin=thin) # samples
* Generating 2000 PMMH samples with 100 particles
  |--------------------------------------------------| 100%
  |**************************************************| 20000 iterations in 2976.07 s

Some summary statistics

summ_pmmh = biips_summary(out_pmmh, probs=c(.025, .975))

Compute kernel density estimates

kde_pmmh = biips_density(out_pmmh)

param_true = log(c_true)
param_lab = expression(log(c[1]), log(c[2]), log(c[3]))

Posterior mean and credible interval of the parameters

for (k in 1:length(param_names)) {
  summ_param = summ_pmmh[[param_names[k]]]
  cat('Posterior mean of ',  param_names[k], ': ', summ_param$mean, '\n', sep='');
  cat('95% credible interval of ',  param_names[k], ': [', summ_param$quant[[1]], ', ', summ_param$quant[[2]],']\n', sep='')
}
Posterior mean of logc[1]: -0.9743927
95% credible interval of logc[1]: [-1.026561, -0.9377507]
Posterior mean of logc[2]: -5.510532
95% credible interval of logc[2]: [-5.595559, -5.451682]
Posterior mean of logc[3]: -0.6937031
95% credible interval of logc[3]: [-0.744729, -0.6284086]

Trace of MCMC samples of the parameters

for (k in 1:length(param_names)) {
  samples_param = out_pmmh[[param_names[k]]]
  plot(samples_param[1,], type='l', col='blue', lwd=1,
       xlab='Iteration', ylab=param_lab[k],
       main=param_lab[k])
  if (sample_data)
    points(0, param_true[k], col='green', pch=8, lwd=2)
}
PMMH: Trace samples parameter

PMMH: Trace samples parameter

Histogram and KDE estimate of the posterior of the parameters

for (k in 1:length(param_names)) {
  samples_param = out_pmmh[[param_names[k]]]
  hist(samples_param, breaks=15, col='blue', border='white',
       xlab=param_lab[k], ylab='Number of samples',
       main=param_lab[k])
  if (sample_data)
    points(param_true[k], 0, col='green', pch=8, lwd=2)
}
PMMH: Histogram posterior parameter

PMMH: Histogram posterior parameter 

for (k in 1:length(param_names)) {
  kde_param = kde_pmmh[[param_names[k]]]
  plot(kde_param, col='blue', lwd=2,
       xlab=param_lab[k], ylab='Posterior density',
       main=param_lab[k])
  if (sample_data)
    points(param_true[k], 0, col='green', pch=8, lwd=2)
}
PMMH: KDE estimate posterior parameter

PMMH: KDE estimate posterior parameter

Posterior mean and quantiles for x

x_pmmh_mean = summ_pmmh$x$mean
x_pmmh_quant = summ_pmmh$x$quant

xx = c(t_vec, rev(t_vec))
yy = c(x_pmmh_quant[[1]][1,], rev(x_pmmh_quant[[2]][1,]))
plot(xx, yy, type='n', xlab='Time', ylab='Number of individuals',
     ylim=c(0, 750))
polygon(xx, yy, col=light_blue, border=NA)

yy = c(x_pmmh_quant[[1]][2,], rev(x_pmmh_quant[[2]][2,]))
polygon(xx, yy, col=light_red, border=NA)

lines(t_vec, x_pmmh_mean[1,], col='blue', lwd=3)
lines(t_vec, data$x_true[1,], col=dark_blue, lwd=2, lty=2)

lines(t_vec, x_pmmh_mean[2,], col='red', lwd=3)
lines(t_vec, data$x_true[2,], col=dark_red, lwd=2, lty=2)

legend('topright', leg=c('95% credible interval (prey)',
                         'PMMH mean estimate (prey)',
                         'True number of preys',
                         '95% credible interval (predator)',
                         'PMMH mean estimate (predator)',
                         'True number of predators'),
       col=c(light_blue,'blue',dark_blue,light_red,'red',dark_red),
       lwd=c(NA,3,2,NA,3,2),
       pch=c(15,NA,NA,15,NA,NA),
       lty=c(NA,1,2,NA,1,2),
       pt.cex=c(2,1,1,2,1,1), bty='n')
SMC: Posterior mean and quantiles

SMC: Posterior mean and quantiles