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

The continuous-time Lotka-Volterra Markov jump process describes the evolution of two species \(X_{1}(t)\) (prey) and \(X_{2}(t)\) (predator) at time \(t\). Let \(dt\) be an infinitesimal interval. The process evolves as

\[\Pr(X_1(t+dt)=x_1(t)+1,X_2(t+dt)=x_2(t)|x_1(t),x_2(t))=c_1x_1(t)dt+o(dt)\]

\[\Pr(X_1(t+dt)=x_1(t)-1,X_2(t+dt)=x_2(t)+1|x_1(t),x_2(t))=c_2x_1(t)x_2(t)dt+o(dt)\]

\[\Pr(X_1(t+dt)=x_1(t),X_2(t+dt)=x_2(t)-1|x_1(t),x_2(t))=c_3 x_2(t)dt+o(dt)\]

where \(c_1=0.5\), \(c_2=0.0025\) and \(c_3=0.3\).

Forward simulation can be done using the Gillespie algorithm. We additionally assume that we observe at some time \(t=1,2,\ldots,t_{\max}\) the number of preys with some noise

\[ Y(t)=X_1(t) + \epsilon(t), ~~\epsilon(t)\sim\mathcal N(0,\sigma^2) \]

Statistical model in BUGS language

Content of the file 'stoch_kinetic_gill.bug':

model_file = 'stoch_kinetic_gill.bug' # BUGS model filename
cat(readLines(model_file), sep = "\n")
# Stochastic kinetic predator-prey model
# cf Boys, Wilkinson and Kirkwood
# Bayesian inference for a discretely observed stochastic kinetic model

var x_true[2,t_max], x[2,t_max], y[t_max], c[3]

data
{
  x_true[,1] ~ LV(x_init,c[1],c[2],c[3],1)
  y[1] ~ dnorm(x_true[1,1], 1/sigma^2) 
  for (t in 2:t_max)
  {      
    x_true[1:2, t] ~ LV(x_true[1:2,t-1],c[1],c[2],c[3],1)   
    y[t] ~ dnorm(x_true[1,t], 1/sigma^2)  
  }
}

model
{
  x[,1] ~ LV(x_init,c[1],c[2],c[3],1)
  y[1] ~ dnorm(x[1,1], 1/sigma^2) 
  for (t in 2:t_max)
  {    
    x[,t] ~ LV(x[,t-1],c[1],c[2],c[3],1) 
    y[t] ~ dnorm(x[1,t], 1/sigma^2) 
  }
}

User-defined R functions

lotka_volterra_gillespie <- function(x, c1, c2, c3, dt) {
  # Simulation from a Lotka-Volterra model with the Gillepsie algorithm
  # x1 is the number of prey
  # x2 is the number of predator
  # R1: (x1,x2) -> (x1+1,x2)      At rate c1x1
  # R2: (x1,x2) -> (x1-1,x2+1)    At rate c2x1x2
  # R3: (x1,x2) -> (x1,x2-1)      At rate c3xx2
  z = matrix(c(1, -1, 0, 0, 1, -1), nrow=2, byrow=TRUE)

  t = 0
  while (TRUE) {
    rate = c(c1*x[1], c2*x[1]*x[2], c3*x[2])
    sum_rate = sum(rate);
    t = t - log(runif(1))/sum_rate # Sample next event from an exponential distribution
    ind = which((sum_rate*runif(1))<=cumsum(rate))[1] # Sample the type of event
    if (t>dt)
      break
    x = x + z[,ind]
  }
  return(x)
}

lotka_volterra_dim <- function(x_dim, c1_dim, c2_dim, c3_dim, dt_dim) { c(2,1) }
  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)

Add new sampler to Biips

Add the user-defined function ‘LV’ to simulate from the Lotka-Volterra model

fun_bugs = 'LV'; fun_nb_inputs = 5
fun_dim = lotka_volterra_dim; fun_sample = lotka_volterra_gillespie
biips_add_distribution(fun_bugs, fun_nb_inputs, fun_dim, fun_sample)
* Added distribution LV

Load model and data

Model parameters

t_max = 40
x_init = c(100, 100)
c = c(.5, .0025, .3)
sigma = 10
data = list(t_max=t_max, c=c, x_init=x_init, sigma=sigma)

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_gill.bug
* Compiling data graph
  Declaring variables
  Resolving undeclared variables
  Allocating nodes
  Graph size: 131
  Sampling data
  Reading data back into data table
* Compiling model graph
  Declaring variables
  Resolving undeclared variables
  Allocating nodes
  Graph size: 132
data = model$data()

Plot data

plot(1:t_max, data$x_true[1,],
     type='l', col='blue', lwd=2,
     xlab='Time',
     ylab='Number of individuals',
     ylim=c(0, 450))
lines(1:t_max, data$x_true[2,], col='red', lwd=2)
points(1:t_max, 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 Sequential Monte Carlo algorithm

Run SMC

n_part = 10000 # Number of particles
variables = c('x') # Variables to be monitored
out_smc = biips_smc_samples(model, variables, n_part, type='fs')
* Assigning node samplers
* Running SMC forward sampler with 10000 particles
  |--------------------------------------------------| 100%
  |**************************************************| 40 iterations in 1310.88 s

Summary statistics

summ_smc = biips_summary(out_smc, probs=c(.025, .975))

Smoothing ESS

plot(out_smc$x$s$ess[1,], type='l', log='y', col='blue', lwd=2,
     xlab='Time', ylab='SESS',
     ylim=c(10, n_part))
lines(1:t_max, rep(30,t_max), lwd=2, lty=2)
SMC: SESS

SMC: SESS

Posterior mean and quantiles for x

x_smc_mean = summ_smc$x$s$mean
x_smc_quant = summ_smc$x$s$quant

xx = c(1:t_max, t_max:1)
yy = c(x_smc_quant[[1]][1,], rev(x_smc_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_smc_quant[[1]][2,], rev(x_smc_quant[[2]][2,]))
polygon(xx, yy, col=light_red, border=NA)

lines(1:t_max, x_smc_mean[1,], col='blue', lwd=3)
lines(1:t_max, data$x_true[1,], col=dark_blue, lwd=2, lty=2)

lines(1:t_max, x_smc_mean[2,], col='red', lwd=3)
lines(1:t_max, data$x_true[2,], col=dark_red, lwd=2, lty=2)

legend('topright', leg=c('95% credible interval (prey)',
                         'SMC mean estimate (prey)',
                         'True number of preys',
                         '95% credible interval (predator)',
                         'SMC 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

Marginal filtering and smoothing density

kde_smc = biips_density(out_smc)
time_index = c(5, 10, 15)
par(mfrow=c(2,2))
for (k in 1:length(time_index)) {
  tk = time_index[k]
  plot(kde_smc$x[[2,tk]], col=c('blue', 'red'), lwd=2,
       xlab=paste('x[2,', tk,']', sep=''), ylab='Posterior density',
       main=paste('Predator at t=', tk, sep=''))
  points(data$x_true[tk], 0, col='green', pch=8, lwd=2)
}
plot(0, type='n', bty='n', xaxt='n', yaxt='n', xlab="", ylab="")
legend('center', leg=c('Filtering density', 'Smoothing density', 'True value'),
       col=c('blue', 'red', 'green'), pch=c(NA,NA,8), lty=c(1,1,NA), lwd=2,
       bty='n')
par(mfrow=c(1,1))
SMC: Marginal posteriors

SMC: Marginal posteriors

Marginal filtering and smoothing probability mass

table_smc = biips_table(out_smc)
time_index = c(5, 10, 15)
par(mfrow=c(2,2))
for (k in 1:length(time_index)) {
  tk = time_index[k]
  plot(table_smc$x[[2,tk]], col=c('blue', 'red'), lwd=2,
       xlab=paste('x[2,', tk,']', sep=''), ylab='Posterior probability mass',
       main=paste('Predator at t=', tk, sep=''))
  points(data$x_true[tk], 0, col='green', pch=8, lwd=2)
}
plot(0, type='n', bty='n', xaxt='n', yaxt='n', xlab="", ylab="")
legend('center', leg=c('Filtering probability mass', 'Smoothing probability mass', 'True value'),
       col=c('blue', 'red', 'green'), pch=c(NA,NA,8), lty=c(1,1,NA), lwd=2,
       bty='n')
par(mfrow=c(1,1))
SMC: Marginal posteriors

SMC: Marginal posteriors