Example simulated time courses from a stochastic Lotka–Volterra model

Description

Collection of simulated time courses from a stochastic Lotka–Volterra model. LVperfect is direct output from a Gillespie simulation. LVprey is the prey component. LVnoise10 has Gaussian noise with standard deviation 10 added. LVnoise30 has Gaussian noise with standard deviation 30 added. LVpreyNoise10 is the prey component with 10 SD noise added. LVnoise3010 has Gaussian noise added. The noise added to the prey component has standard deviation 30 and the noise added to the predator component has standard deviation 10. LVnoise10scale10 has Gaussian noise with standard deviation 10 added, and is then rescaled by a factor of 10 to mimic a scenario of an uncalibrated measurement scale. LVirregular is direct output from a Gillespie simulator, but on an irregular time grid. LVirregularNoise10 is output on an irregular time grid with Gaussian noise of standard deviation 10 added.

Usage

data(LVdata)

Format

All datasets beginning LVirregular are R matrices such as output by simTimes, and the rest are R ts objects such as output by simTs.

Create a function for advancing the state of an SPN by using a simple Euler-Maruyama integration method for the approximating CLE

Description

This function creates a function for advancing the state of an SPN model using a simple Euler-Maruyama integration method for the approximating chemical Langevin equation (CLE). The resulting function (closure) can be used in conjunction with other functions (such as simTs) for simulating realisations of SPN models.

Usage

StepCLE(N,dt=0.01)

Arguments

Value

An R function which can be used to advance the state of the SPN model N by using an Euler-Maruyama method on the approximating CLE with step size dt. The function closure has interface function(x0,t0,deltat,...), where x0 and t0 represent the initial state and time, and deltat represents the amount of time by which the process should be advanced. The function closure returns a vector representing the simulated state of the system at the new time.

See Also

StepGillespie, StepEulerSPN, StepSDE, simTs, simSample

Examples

# load the LV model
data(spnModels)
# create a stepping function
stepLV = StepCLE(LV)
# step the function
print(stepLV(c(x1=50,x2=100),0,1))
# integrate the process and plot it
out = simTs(c(x1=50,x2=100),0,20,0.1,stepLV)
plot(out)
plot(out,plot.type="single",lty=1:2)

Create a function for advancing the state of an SPN by using a simple Euler-Maruyama discretisation of the CLE on a 1D regular grid

Description

This function creates a function for advancing the state of an SPN model using a simple Euler-Maruyama discretisation of the CLE on a 1D regular grid. The resulting function (closure) can be used in conjunction with other functions (such as simTs1D) for simulating realisations of SPN models in space and time.

Usage

StepCLE1D(N,d,dt=0.01)

Arguments

Value

An R function which can be used to advance the state of the SPN model N by using a simple Euler-Maruyama algorithm. The function closure has interface function(x0,t0,deltat,...), where x0 is a matrix with rows corresponding to species and columns corresponding to voxels, representing the initial condition, t0 represent the initial state and time, and deltat represents the amount of time by which the process should be advanced. The function closure returns a matrix representing the simulated state of the system at the new time.

See Also

StepGillespie1D,StepCLE, simTs1D, StepCLE2D

Examples

N=200
T=40
data(spnModels)
x0=matrix(0,nrow=2,ncol=N)
rownames(x0)=c("x1","x2")
x0[,round(N/2)]=LV$M
stepLV1D = StepCLE1D(LV,c(0.6,0.6),dt=0.05)
xx = simTs1D(x0,0,T,0.2,stepLV1D)
op=par(mfrow=c(1,2))
image(xx[1,,],main="Prey",xlab="Space",ylab="Time")
image(xx[2,,],main="Predator",xlab="Space",ylab="Time")
par(op)

Create a function for advancing the state of an SPN by using a simple Euler-Maruyama discretisation of the CLE on a 2D regular grid

Description

This function creates a function for advancing the state of an SPN model using a simple Euler-Maruyama discretisation of the CLE on a 2D regular grid. The resulting function (closure) can be used in conjunction with other functions (such as simTs2D) for simulating realisations of SPN models in space and time.

Usage

StepCLE2D(N,d,dt=0.01)

Arguments

Value

An R function which can be used to advance the state of the SPN model N by using a simple Euler-Maruyama algorithm. The function closure has interface function(x0,t0,deltat,...), where x0 is a 3D array with rows corresponding to species and columns corresponding to voxels, representing the initial condition (with dimensions species, x, and y), t0 represent the initial state and time, and deltat represents the amount of time by which the process should be advanced. The function closure returns a matrix representing the simulated state of the system at the new time.

See Also

StepGillespie2D,StepCLE, simTs1D, StepCLE1D

Examples

m=150
n=100
T=15
data(spnModels)
x0=array(0,c(2,m,n))
dimnames(x0)[[1]]=c("x1","x2")
x0[,round(m/2),round(n/2)]=LV$M
stepLV2D = StepCLE2D(LV,c(0.6,0.6),dt=0.05)
xx = simTs2D(x0,0,T,0.5,stepLV2D,verb=TRUE)
N = dim(xx)[4]
op=par(mfrow=c(1,2))
image(xx[1,,,N],main="Prey",xlab="Space",ylab="Time")
image(xx[2,,,N],main="Predator",xlab="Space",ylab="Time")
par(op)

Create a function for advancing the state of an ODE model by using a simple Euler integration method

Description

This function creates a function for advancing the state of an ODE model using a simple Euler integration method. The resulting function (closure) can be used in conjunction with other functions (such as simTs) for simulating realisations of ODE models. This function is intended to be pedagogic. See StepODE for a more accurate integration function.

Usage

StepEuler(RHSfun,dt=0.01)

Arguments

Value

An R function which can be used to advance the state of the ODE model RHSfun by using an Euler method with step size dt. The function closure has interface function(x0,t0,deltat,...), where t0 and x0 represent the initial time and state, and deltat represents the amount of time by which the process should be advanced. The function closure returns a vector representing the simulated state of the system at the new time.

See Also

StepEulerSPN, StepODE, simTs, simSample

Examples

# Build a RHS for the Lotka-Volterra system
LVrhs <- function(x,t,th=c(c1=1,c2=0.005,c3=0.6))
{
        with(as.list(c(x,th)),{
                c( c1*x1 - c2*x1*x2 ,
                      c2*x1*x2 - c3*x2 )
        })
}
# create a stepping function
stepLV = StepEuler(LVrhs)
# step the function
print(stepLV(c(x1=50,x2=100),0,1))
# integrate the process and plot it
out = simTs(c(x1=50,x2=100),0,20,0.1,stepLV)
plot(out,plot.type="single",lty=1:2)

Create a function for advancing the state of an SPN by using a simple continuous deterministic Euler integration method

Description

This function creates a function for advancing the state of an SPN model using a simple continuous deterministic Euler integration method. The resulting function (closure) can be used in conjunction with other functions (such as simTs) for simulating realisations of SPN models.

Usage

StepEulerSPN(N,dt=0.01)

Arguments

Value

An R function which can be used to advance the state of the SPN model N by using an Euler method with step size dt. The function closure has interface function(x0,t0,deltat,...), where x0 and t0 represent the initial state and time, and deltat represents the amount of time by which the process should be advanced. The function closure returns a vector representing the simulated state of the system at the new time.

See Also

StepGillespie, StepODE, StepCLE, simpleEuler, simTs, simSample

Examples

# load the LV model
data(spnModels)
# create a stepping function
stepLV = StepEulerSPN(LV)
# step the function
print(stepLV(c(x1=50,x2=100),0,1))
# integrate the process and plot it
out = simTs(c(x1=50,x2=100),0,100,0.1,stepLV)
plot(out)
plot(out,plot.type="single",lty=1:2)

Create a function for advancing the state of an SPN by using Gillespie's first reaction method

Description

This function creates a function for advancing the state of an SPN model using Gillespie's first reaction method. The resulting function (closure) can be used in conjunction with other functions (such as simTs) for simulating realisations of SPN models.

Usage

StepFRM(N)

Arguments

Value

An R function which can be used to advance the state of the SPN model N by using Gillespie's first reaction method. The function closure has interface function(x0,t0,deltat,...), where x0 and t0 represent the initial state and time, and deltat represents the amount of time by which the process should be advanced. The function closure returns a vector representing the simulated state of the system at the new time.

See Also

StepEulerSPN, StepGillespie, simTs, simSample

Examples

# load the LV model
data(spnModels)
# create a stepping function
stepLV = StepFRM(LV)
# step the function
print(stepLV(c(x1=50,x2=100),0,1))
# simulate a realisation of the process and plot it
out = simTs(c(x1=50,x2=100),0,100,0.1,stepLV)
plot(out,plot.type="single",lty=1:2)

Create a function for advancing the state of an SPN by using the Gillespie algorithm

Description

This function creates a function for advancing the state of an SPN model using the Gillespie algorithm. The resulting function (closure) can be used in conjunction with other functions (such as simTs) for simulating realisations of SPN models.

Usage

StepGillespie(N)

Arguments

Value

An R function which can be used to advance the state of the SPN model N by using the Gillespie algorithm. The function closure has interface function(x0,t0,deltat,...), where x0 and t0 represent the initial state and time, and deltat represents the amount of time by which the process should be advanced. The function closure returns a vector representing the simulated state of the system at the new time.

See Also

StepEulerSPN, StepGillespie1D, simTs, simTimes, simSample, StepFRM, StepPTS, StepCLE

Examples

# load up the Lotka-Volterra (LV) model
data(spnModels)
LV
# create a stepping function
stepLV = StepGillespie(LV)
# step the function
print(stepLV(c(x1=50,x2=100),0,1))
# simulate a realisation of the process and plot it
out = simTs(c(x1=50,x2=100),0,100,0.1,stepLV)
plot(out)
plot(out,plot.type="single",lty=1:2)
# simulate a realisation using simTimes
times = seq(0,100,by=0.1)
plot(ts(simTimes(c(x1=50,x2=100),0,times,stepLV),start=0,deltat=0.1),plot.type="single",lty=1:2)
# simulate a realisation at irregular times
times = c(0,10,20,50,100)
out2 = simTimes(c(x1=50,x2=100),0,times,stepLV)
print(out2)

Create a function for advancing the state of an SPN by using the Gillespie algorithm on a 1D regular grid

Description

This function creates a function for advancing the state of an SPN model using the Gillespie algorithm. The resulting function (closure) can be used in conjunction with other functions (such as simTs1D) for simulating realisations of SPN models in space and time.

Usage

StepGillespie1D(N,d)

Arguments

Value

An R function which can be used to advance the state of the SPN model N by using the Gillespie algorithm. The function closure has interface function(x0,t0,deltat,...), where x0 is a matrix with rows corresponding to species and columns corresponding to voxels, representing the initial condition, t0 represent the initial state and time, and deltat represents the amount of time by which the process should be advanced. The function closure returns a matrix representing the simulated state of the system at the new time.

See Also

StepGillespie, simTs1D, StepGillespie2D

Examples

data(spnModels)
N=20; T=30
x0=matrix(0,nrow=2,ncol=N)
rownames(x0)=c("x1","x2")
x0[,round(N/2)]=LV$M
stepLV1D = StepGillespie1D(LV,c(0.6,0.6))
xx = simTs1D(x0,0,T,0.2,stepLV1D,verb=TRUE)
op=par(mfrow=c(1,2))
image(xx[1,,],main="Prey",xlab="Space",ylab="Time")
image(xx[2,,],main="Predator",xlab="Space",ylab="Time")
par(op)

Create a function for advancing the state of an SPN by using the Gillespie algorithm on a 2D regular grid

Description

This function creates a function for advancing the state of an SPN model using the Gillespie algorithm. The resulting function (closure) can be used in conjunction with other functions (such as simTs2D) for simulating realisations of SPN models in space and time.

Usage

StepGillespie2D(N,d)

Arguments

Value

An R function which can be used to advance the state of the SPN model N by using the Gillespie algorithm. The function closure has interface function(x0,t0,deltat,...), where x0 is a 3d array with dimensions corresponding to species followed by two spatial dimensions, representing the initial condition, t0 represent the initial state and time, and deltat represents the amount of time by which the process should be advanced. The function closure returns an array representing the simulated state of the system at the new time.

See Also

StepGillespie, simTs2D, StepGillespie1D

Examples

data(spnModels)
m=20; n=30; T=10
x0=array(0,c(2,m,n))
dimnames(x0)[[1]]=c("x1","x2")
x0[,round(m/2),round(n/2)]=LV$M
stepLV2D = StepGillespie2D(LV,c(0.6,0.6))
xx = simTs2D(x0,0,T,0.2,stepLV2D,verb=TRUE)
N = dim(xx)[4]
op=par(mfrow=c(1,2))
image(xx[1,,,N],main="Prey",xlab="Space",ylab="Time")
image(xx[2,,,N],main="Predator",xlab="Space",ylab="Time")
par(op)

Create a function for advancing the state of an ODE model by using the deSolve package

Description

This function creates a function for advancing the state of an ODE model using an integration method from the deSolve package. The resulting function (closure) can be used in conjunction with other functions (such as simTs) for simulating realisations of ODE models. This function is used similarly to StepEuler, but StepODE should be more accurate and efficient.

Usage

StepODE(RHSfun)

Arguments

Value

An R function which can be used to advance the state of the ODE model RHSfun by using an efficient ODE solver. The function closure has interface function(x0,t0,deltat,parms,...), where t0 and x0 represent the initial time and state, and deltat represents the amount of time by which the process should be advanced. The function closure returns a vector representing the simulated state of the system at the new time.

See Also

StepEulerSPN, StepEuler, simTs, ode

Examples

# Build a RHS for the Lotka-Volterra system
LVrhs <- function(x,t,parms)
{
        with(as.list(c(x,parms)),{
                c( c1*x1 - c2*x1*x2 ,
                      c2*x1*x2 - c3*x2 )
        })
}
# create a stepping function
stepLV = StepODE(LVrhs)
# step the function
print(stepLV(c(x1=50,x2=100),0,1,parms=c(c1=1,c2=0.005,c3=0.6)))
# integrate the process and plot it
out = simTs(c(x1=50,x2=100),0,50,0.1,stepLV,parms=c(c1=1,c2=0.005,c3=0.6))
plot(out,plot.type="single",lty=1:2)

Create a function for advancing the state of an SPN by using a simple approximate Poisson time stepping method

Description

This function creates a function for advancing the state of an SPN model using a simple approximate Poisson time stepping method. The resulting function (closure) can be used in conjunction with other functions (such as simTs) for simulating realisations of SPN models.

Usage

StepPTS(N,dt=0.01)

Arguments

Value

An R function which can be used to advance the state of the SPN model N by using a Poisson time stepping method with step size dt. The function closure has interface function(x0,t0,deltat,...), where x0 and t0 represent the initial state and time, and deltat represents the amount of time by which the process should be advanced. The function closure returns a vector representing the simulated state of the system at the new time.

See Also

StepGillespie, StepCLE, simTs, simSample

Examples

# load up the LV model
data(spnModels)
# create a stepping function
stepLV=StepPTS(LV)
# step the function
print(stepLV(c(x1=50,x2=100),0,1))
# integrate the process and plot it
out = simTs(c(x1=50,x2=100),0,20,0.1,stepLV)
plot(out)
plot(out,plot.type="single",lty=1:2)

Create a function for advancing the state of an SDE model by using a simple Euler-Maruyama integration method

Description

This function creates a function for advancing the state of an SDE model using a simple Euler-Maruyama integration method. The resulting function (closure) can be used in conjunction with other functions (such as simTs) for simulating realisations of SDE models.

Usage

StepSDE(drift,diffusion,dt=0.01)

Arguments

Value

An R function which can be used to advance the state of the SDE model with given drift vector and diffusion matrix, by using an Euler-Maruyama method with step size dt. The function closure has interface function(x0,t0,deltat,...), where x0 and t0 represent the initial state and time, and deltat represents the amount of time by which the process should be advanced. The function closure returns a vector representing the simulated state of the system at the new time.

See Also

StepEuler, StepCLE, simTs, simSample

Examples

# Immigration-death diffusion approx with death rate a CIR process
myDrift <- function(x,t,th=c(lambda=1,alpha=1,mu=0.1,sigma=0.1))
     {
             with(as.list(c(x,th)),{
                     c( lambda - x*y ,
                           alpha*(mu-y) )
             })
     }
myDiffusion <- function(x,t,th=c(lambda=1,alpha=1,mu=0.1,sigma=0.1))
     {
             with(as.list(c(x,th)),{
                     matrix(c( sqrt(lambda + x*y) , 0,
                           0, sigma*sqrt(y) ),ncol=2,nrow=2,byrow=TRUE)
             })
     }
# create a stepping function
stepProc = StepSDE(myDrift,myDiffusion)
# integrate the process and plot it
out = simTs(c(x=5,y=0.1),0,20,0.1,stepProc)
plot(out)

Run a set of simulations initialised with parameters sampled from a given prior distribution, and compute statistics required for an ABC analaysis

Description

Run a set of simulations initialised with parameters sampled from a given prior distribution, and compute statistics required for an ABC analaysis. Typically used to calculate "distances" of simulated synthetic data from observed data.

Usage

abcRun(n, rprior, rdist)

Arguments

Value

A list with elements 'param' and 'dist'. These will be returned as matrices or vectors depending on whether the parameters and distances are scalars or vectors.

See Also

pfMLLik, StepGillespie, abcSmc, simTs, stepLVc

Examples

data(LVdata)
rprior <- function() { exp(c(runif(1, -3, 3),runif(1,-8,-2),runif(1,-4,2))) }
rmodel <- function(th) { simTs(c(50,100), 0, 30, 2, stepLVc, th) }
sumStats <- identity
ssd = sumStats(LVperfect)
distance <- function(s) {
    diff = s - ssd
    sqrt(sum(diff*diff))
}
rdist <- function(th) { distance(sumStats(rmodel(th))) }
out = abcRun(10000, rprior, rdist)
q=quantile(out$dist, c(0.01, 0.05, 0.1))
print(q)
accepted = out$param[out$dist < q[1],]
print(summary(accepted))
print(summary(log(accepted)))

Run an ABC-SMC algorithm for infering the parameters of a forward model

Description

Run an ABC-SMC algorithm for infering the parameters of a forward model. This sequential Monte Carlo algorithm often performs better than simple rejection-ABC in practice.

Usage

abcSmc(N, rprior, dprior, rdist, rperturb, dperturb, factor=10,
                   steps=15, verb=FALSE)

Arguments

Value

A matrix (or vector) with rows (or elements) representing samples from the approximate posterior distribution.

See Also

pfMLLik, StepGillespie, abcRun, simTs, stepLVc

Examples

data(LVdata)
rprior <- function() { c(runif(1, -3, 3), runif(1, -8, -2), runif(1, -4, 2)) }
dprior <- function(x, ...) { dunif(x[1], -3, 3, ...) + 
                dunif(x[2], -8, -2, ...) + dunif(x[3], -4, 2, ...) }
rmodel <- function(th) { simTs(c(50,100), 0, 30, 2, stepLVc, exp(th)) }
rperturb <- function(th){th + rnorm(3, 0, 0.5)}
dperturb <- function(thNew, thOld, ...){sum(dnorm(thNew, thOld, 0.5, ...))}
sumStats <- identity
ssd = sumStats(LVperfect)
distance <- function(s) {
    diff = s - ssd
    sqrt(sum(diff*diff))
}
rdist <- function(th) { distance(sumStats(rmodel(th))) }
out = abcSmc(5000, rprior, dprior, rdist, rperturb,
             dperturb, verb=TRUE, steps=6, factor=5)
print(summary(out))

Convert a time series object to a timed data matrix

Description

This function converts a time series object to a timed data matrix, similar to that produced by simTimes. The main purpose is for passing data to the function pfMLLik, which expects data encoded in this format.

Usage

as.timedData(timeseries)

Arguments

Value

An R matrix object with row names corresponding to observation times, similar to that produced by simTimes.

See Also

simTs, ts, simTimes, pfMLLik

Examples

truth=simTs(c(x1=50,x2=100),0,20,2,stepLVc)
simData=truth+rnorm(prod(dim(truth)),0,5)
timedData=as.timedData(simData)
print(timedData)

Discretise output from a discrete event simulation algorithm

Description

This function discretise output from a discrete event simulation algorithm such as gillespie onto a regular time grid, and returns the results as an R ts object.

Usage

discretise(out, dt=1, start=0)

Arguments

Value

An R ts object containing the discretised output.

See Also

simpleEuler, rdiff, gillespie, gillespied, ts

Examples

# load LV model
data(spnModels)
# simulate a realisation of the process and plot it
out = gillespie(LV,10000)
op=par(mfrow=c(2,2))
plot(stepfun(out$t,out$x[,1]),pch="")
plot(stepfun(out$t,out$x[,2]),pch="")
plot(out$x,type="l")

# use the "discretise" function to map it to an R "ts" object
plot(discretise(out,dt=0.01),plot.type="single",lty=c(1,2))
par(op)

Simulate a sample path from a stochastic kinetic model described by a stochastic Petri net

Description

This function simulates a single realisation from a discrete stochastic kinetic model described by a stochastic Petri net (SPN).

Usage

gillespie(N, n, ...)

Arguments

Value

A list with first component t, a vector of length n containing event times and second component x, a matrix with n+1 rows containing the state of the system. The ith row of x contains the state of the system prior to the ith event.

See Also

simpleEuler, rdiff, discretise, gillespied, StepGillespie

Examples

# load the LV model
data(spnModels)
# simulate a realisation of the process and plot it
out = gillespie(LV,10000)
op = par(mfrow=c(2,2))
plot(stepfun(out$t,out$x[,1]),pch="")
plot(stepfun(out$t,out$x[,2]),pch="")
plot(out$x,type="l")

# use the "discretise" function to map it to an R "ts" object
plot(discretise(out,dt=0.01),plot.type="single",lty=c(1,2))
par(op)

Simulate a sample path from a stochastic kinetic model described by a stochastic Petri net

Description

This function simulates a single realisation from a discrete stochastic kinetic model described by a stochastic Petri net and discretises the output onto a regular time grid.

Usage

gillespied(N, T=100, dt=1, ...)

Arguments

Value

An R ts object containing the simulated realisation of the process.

See Also

simpleEuler, rdiff, discretise, gillespie, StepGillespie

Examples

# load LV model
data(spnModels)
# simulate and plot a realisation
plot(gillespied(LV,T=100,dt=0.01))

Simulate a sample path from the homogeneous immigration-death process

Description

This function simulates a single realisation from a time-homogeneous immigration-death process.

Usage

imdeath(n=20,x0=0,lambda=1,mu=0.1)

Arguments

Value

An R stepfun object containing the sampled path of the process.

See Also

rcfmc, rdiff, stepfun, gillespie

Examples

plot(imdeath(50))

Summarise and plot tabular MCMC output

Description

This function summarises and plots tabular MCMC output such as that generated by the function normgibbs.

Usage

mcmcSummary(mat, rows = 4, lag.max=100, bins=30, show = TRUE, plot = TRUE, truth = NULL) 

Arguments

Value

An R summary object.

See Also

normgibbs, summary, acf

Examples

out=normgibbs(N=1000,n=15,a=3,b=11,cc=10,d=1/100,xbar=25,ssquared=20)
names(out)=c("mu","tau")
mcmcSummary(out,rows=2,bins=10,truth=c(25,1/20))

Run a simple Metropolis sampler with standard normal target and uniform innovations

Description

This function runs a simple Metropolis sampler with standard normal target distribution and uniform innovations.

Usage

metrop(n, alpha) 

Arguments

Value

An R vector containing the output of the sampler.

See Also

normgibbs

Examples

normvec=metrop(1000,1)
op=par(mfrow=c(2,1))
plot(ts(normvec))
hist(normvec,20)
par(op)

Run a Metropolis-Hastings MCMC algorithm for the parameters of a Bayesian posterior distribution

Description

Run a Metropolis-Hastings MCMC algorithm for the parameters of a Bayesian posterior distribution. Note that the algorithm carries over the old likelihood from the previous iteration, making it suitable for problems with expensive likelihoods, and also for "exact approximate" pseudo-marginal or particle marginal MH algorithms.

Usage

metropolisHastings(init, logLik, rprop, dprop=function(new, old, ...){1},
                   dprior=function(x, ...){1}, iters=10000, thin=10,
                   verb=TRUE, debug=FALSE)

Arguments

Value

A matrix with rows representing samples from the posterior distribution.

See Also

pfMLLik, StepGillespie, abcRun, simTs, stepLVc, metrop

Examples

## First simulate some synthetic data
data = rnorm(250,5,2)
## Now use MH to recover the parameters
llik = function(x) { sum(dnorm(data,x[1],x[2],log=TRUE)) }
prop = function(x) { rnorm(2,x,0.1) }
prior = function(x, log=TRUE) {
    l = dnorm(x[1],0,100,log=TRUE) + dgamma(x[2],1,0.0001,log=TRUE)
    if (log) l else exp(l)
}
out = metropolisHastings(c(mu=1,sig=1), llik, prop,
                         dprior=prior, verb=FALSE)
out = out[1000:10000,]
mcmcSummary(out, truth=c(5,2), rows=2, plot=FALSE)

Simple example data frame

Description

Trivial example of a very small data frame. Used as part of the R tutorial.

Usage

data(mytable)

Format

A very small example data frame.

A simple Gibbs sampler for Bayesian inference for the mean and precision of a normal random sample

Description

This function runs a simple Gibbs sampler for the Bayesian posterior distribution of the mean and precision given a normal random sample.

Usage

normgibbs(N, n, a, b, cc, d, xbar, ssquared)

Arguments

Value

An R matrix object containing the samples of the Gibbs sampler.

See Also

rcfmc, metrop, mcmcSummary

Examples

postmat=normgibbs(N=1100,n=15,a=3,b=11,cc=10,d=1/100,xbar=25,ssquared=20)
postmat=postmat[101:1100,]
op=par(mfrow=c(3,3))
plot(postmat)
plot(postmat,type="l")
plot.new()
plot(ts(postmat[,1]))
plot(ts(postmat[,2]))
plot(ts(sqrt(1/postmat[,2])))
hist(postmat[,1],30)
hist(postmat[,2],30)
hist(sqrt(1/postmat[,2]),30)
par(op)

Create a function for computing the log of an unbiased estimate of marginal likelihood of a time course data set

Description

Create a function for computing the log of an unbiased estimate of marginal likelihood of a time course data set using a simple bootstrap particle filter. This version uses the "log-sum-exp trick" for avoiding numerical underflow of weights. See pfMLLik1 for a version which doesn't.

Usage

pfMLLik(n,simx0,t0,stepFun,dataLik,data)

Arguments

Value

An R function with interface (...) which evaluates to the log of the particle filter's unbiased estimate of the marginal likelihood of the data.

See Also

pfMLLik1, StepGillespie, as.timedData, simTimes, stepLVc

Examples

noiseSD=5
# first simulate some data
truth=simTs(c(x1=50,x2=100),0,20,2,stepLVc)
data=truth+rnorm(prod(dim(truth)),0,noiseSD)
data=as.timedData(data)
# measurement error model
dataLik <- function(x,t,y,log=TRUE,...)
{
    ll=sum(dnorm(y,x,noiseSD,log=TRUE))
    if (log)
        return(ll)
    else
        return(exp(ll))
}
# now define a sampler for the prior on the initial state
simx0 <- function(N,t0,...)
{
    mat=cbind(rpois(N,50),rpois(N,100))
    colnames(mat)=c("x1","x2")
    mat
}
mLLik=pfMLLik(1000,simx0,0,stepLVc,dataLik,data)
print(mLLik())
print(mLLik(th=c(th1 = 1, th2 = 0.005, th3 = 0.6)))
print(mLLik(th=c(th1 = 1, th2 = 0.005, th3 = 0.5)))

Create a function for computing the log of an unbiased estimate of marginal likelihood of a time course data set

Description

Create a function for computing the log of an unbiased estimate of marginal likelihood of a time course data set using a simple bootstrap particle filter. This version does not use the "log-sum-exp trick" for avoiding numerical underflow. See pfMLLik for a version which does.

Usage

pfMLLik1(n,simx0,t0,stepFun,dataLik,data)

Arguments

Value

An R function with interface (...) which evaluates to the log of the particle filter's unbiased estimate of the marginal likelihood of the data.

See Also

pfMLLik, StepGillespie, as.timedData, simTimes, stepLVc

Examples

noiseSD=5
# first simulate some data
truth=simTs(c(x1=50,x2=100),0,20,2,stepLVc)
data=truth+rnorm(prod(dim(truth)),0,noiseSD)
data=as.timedData(data)
# measurement error model
dataLik <- function(x,t,y,log=TRUE,...)
{
    ll=sum(dnorm(y,x,noiseSD,log=TRUE))
    if (log)
        return(ll)
    else
        return(exp(ll))
}
# now define a sampler for the prior on the initial state
simx0 <- function(N,t0,...)
{
    mat=cbind(rpois(N,50),rpois(N,100))
    colnames(mat)=c("x1","x2")
    mat
}
mLLik=pfMLLik1(1000,simx0,0,stepLVc,dataLik,data)
print(mLLik())
print(mLLik(th=c(th1 = 1, th2 = 0.005, th3 = 0.6)))
print(mLLik(th=c(th1 = 1, th2 = 0.005, th3 = 0.5)))

Simulate a continuous time finite state space Markov chain

Description

This function simulates a single realisation from a continuous time Markov chain having a finite state space based on a given transition rate matrix.

Usage

rcfmc(n,Q,pi0)

Arguments

Value

An R stepfun object containing the sampled path of the process.

See Also

rfmc, stepfun

Examples

plot(rcfmc(20,matrix(c(-0.5,0.5,1,-1),ncol=2,byrow=TRUE),c(1,0)))

Simulate a sample path from a univariate diffusion process

Description

This function simulates a single realisation from a time-homogeneous univariate diffusion process.

Usage

rdiff(afun, bfun, x0 = 0, t = 50, dt = 0.01, ...)

Arguments

Value

An R ts object containing the sampled path of the process.

See Also

rcfmc, ts, StepSDE, simTs

Examples

# simulate a diffusion approximation to an immigration-death process
# infinitesimal mean
afun<-function(x,lambda,mu)
{
        lambda-mu*x
}
# infinitesimal standard deviation
bfun<-function(x,lambda,mu)
{
        sqrt(lambda+mu*x)
}
# plot a sample path
plot(rdiff(afun,bfun,lambda=1,mu=0.1,t=30))

Simulate a finite state space Markov chain

Description

This function simulates a single realisation from a discrete time Markov chain having a finite state space based on a given transition matrix.

Usage

rfmc(n,P,pi0)

Arguments

Value

An R ts object containing the sampled values from the Markov chain.

See Also

rcfmc, ts

Examples

# example for sampling a finite Markov chain
P = matrix(c(0.9,0.1,0.2,0.8),ncol=2,byrow=TRUE)
pi0 = c(0.5,0.5)
samplepath = rfmc(200,P,pi0)
plot(samplepath)
summary(samplepath)
table(samplepath)
table(samplepath)/length(samplepath) # empirical distribution
# now compute the exact stationary distribution...
e = eigen(t(P))$vectors[,1]
e/sum(e)

Simulate a many realisations of a model at a given fixed time in the future given an initial time and state, using a function (closure) for advancing the state of the model

Description

This function simulates many realisations of a model at a given fixed time in the future given an initial time and state, using a function (closure) for advancing the state of the model , such as created by StepGillespie or StepSDE.

Usage

simSample(n=100,x0,t0=0,deltat,stepFun,...)

Arguments

Value

An R matrix whose rows represent the simulated states of the process at time t0+deltat.

See Also

StepSDE, StepGillespie, simTimes, simTs

Examples

out3 = simSample(100,c(x1=50,x2=100),0,20,stepLVc)
hist(out3[,"x2"])

Simulate a model at a specified set of times, using a function (closure) for advancing the state of the model

Description

This function simulates a single realisation from a Markovian model and records the state at a specified set of times using a function (closure) for advancing the state of the model, such as created by StepGillespie or StepEulerSPN.

Usage

simTimes(x0,t0=0,times,stepFun,...)

Arguments

Value

An R matrix where each row represents the state of the process at one of the required times. The row names contain the sampled times.

See Also

StepEulerSPN, StepGillespie, simTs, simSample, as.timedData, pfMLLik

Examples

# load the LV model
data(spnModels)
# create a stepping function
stepLV = StepGillespie(LV)
# simulate a realisation using simTimes
times = seq(0,100,by=0.1)
plot(ts(simTimes(c(x1=50,x2=100),0,times,stepLV),start=0,deltat=0.1),plot.type="single",lty=1:2)
# simulate a realisation at irregular times
times = c(0,10,20,50,100)
out2 = simTimes(c(x1=50,x2=100),0,times,stepLV)
print(out2)

Simulate a model on a regular grid of times, using a function (closure) for advancing the state of the model

Description

This function simulates single realisation of a model on a regular grid of times using a function (closure) for advancing the state of the model, such as created by StepGillespie or StepEulerSPN.

Usage

simTs(x0,t0=0,tt=100,dt=0.1,stepFun,...)

Arguments

Value

An R ts object representing the simulated process.

See Also

StepEulerSPN, StepGillespie, StepSDE, simTimes, simSample, as.timedData

Examples

# load the LV model
data(spnModels)
# create a stepping function
stepLV = StepGillespie(LV)
# simulate a realisation of the process and plot it
out = simTs(c(x1=50,x2=100),0,100,0.1,stepLV)
plot(out)
plot(out,plot.type="single",lty=1:2)

Simulate a model on a regular grid of times, using a function (closure) for advancing the state of the model

Description

This function simulates single realisation of a model on a 1D regular spatial grid and regular grid of times using a function (closure) for advancing the state of the model, such as created by StepGillespie1D.

Usage

simTs1D(x0,t0=0,tt=100,dt=0.1,stepFun,verb=FALSE,...)

Arguments

Value

An R 3d array representing the simulated process. The dimensions are species, space, and time.

See Also

StepGillespie1D, simTs

Examples

data(spnModels)
N=20; T=30
x0=matrix(0,nrow=2,ncol=N)
rownames(x0)=c("x1","x2")
x0[,round(N/2)]=LV$M
stepLV1D = StepGillespie1D(LV,c(0.6,0.6))
xx = simTs1D(x0,0,T,0.2,stepLV1D,verb=TRUE)
op=par(mfrow=c(1,2))
image(xx[1,,],main="Prey",xlab="Space",ylab="Time")
image(xx[2,,],main="Predator",xlab="Space",ylab="Time")
par(op)

Simulate a model on a regular grid of times, using a function (closure) for advancing the state of the model

Description

This function simulates single realisation of a model on a 2D regular spatial grid and regular grid of times using a function (closure) for advancing the state of the model, such as created by StepGillespie2D.

Usage

simTs2D(x0,t0=0,tt=100,dt=0.1,stepFun,verb=FALSE,...)

Arguments

Value

An R 4d array representing the simulated process. The dimensions are species, 2 space, and time.

See Also

StepGillespie2D, simTs1D

Examples

data(spnModels)
m=20; n=30; T=15
x0=array(0,c(2,m,n))
dimnames(x0)[[1]]=c("x1","x2")
x0[,round(m/2),round(n/2)]=LV$M
stepLV2D = StepGillespie2D(LV,c(0.6,0.6))
xx = simTs2D(x0,0,T,0.2,stepLV2D,verb=TRUE)
N = dim(xx)[4]
op=par(mfrow=c(1,2))
image(xx[1,,,N],main="Prey",xlab="Space",ylab="Time")
image(xx[2,,,N],main="Predator",xlab="Space",ylab="Time")
par(op)

Simulate a sample path from an ODE model

Description

This function integrates an Ordinary Differential Equation (ODE) model using a simple first order Euler method. The function is pedagogic and not intended for serious use. See the deSolve package for better, more robust ODE solvers.

Usage

simpleEuler(t=50, dt=0.001, fun, ic, ...)

Arguments

Value

An R ts object containing the sampled path of the model.

See Also

rdiff, ts, StepEuler, simTs

Examples

# simple Lotka-Volterra example
lv <- function(x,t,k=c(k1=1,k2=0.1,k3=0.1))
{
        with(as.list(c(x,k)),{
                c( k1*x1 - k2*x1*x2 ,
                      k2*x1*x2 - k3*x2 )
        })
}
plot(simpleEuler(t=100,fun=lv,ic=c(x1=4,x2=10)),plot.type="single",lty=1:2)

# now an example which instead uses deSolve...
require(deSolve)
times = seq(0,50,by=0.01)
k = c(k1=1,k2=0.1,k3=0.1)
lvlist = function(t,x,k)
        list(lv(x,t,k))
plot(ode(y=c(x1=4,x2=10),times=times,func=lvlist,parms=k))

Stochastic Modelling for Systems Biology

Description

This package contains code and data for modelling and simulation of stochastic kinetic biochemical network models. It contains the code and data associated with the second and third editions of the book Stochastic Modelling for Systems Biology, published by Chapman & Hall/CRC Press.

Author(s)

Maintainer: Darren Wilkinson <darrenjwilkinson@btinternet.com>

References

See https://darrenjw.github.io/work/smfsb/ or https://github.com/darrenjw/smfsb for further details.

Example SPN models

Description

Collection of example stochastic Petri net (SPN) models. Includes LV, a Lotka–Volterra model, ID, an immigration–death process, BD, a birth–death process, SIR, a simple SIR model, SEIR, an SEIR epidemic model, Dimer, a simple dimerisation kinetics model, and MM, a Michaelis–Menten enzyme kinetic model.

Usage

data(spnModels)

Format

Each model is a list, with components Pre, Post, and h. Some models also include an initial state, M. See gillespie and StepGillespie for further details, and examples of use.

A function for advancing the state of a Lotka-Volterra model by using the Gillespie algorithm

Description

A function for advancing the state of a Lotka-Volterra model by calling some C code implementing the Gillespie algorithm. The function can be used in conjunction with other functions (such as simTs) for simulating realisations of Lotka-Volterra models. Should be functionally identical to the function obtained by data(spnModels), stepLV=StepGillespie(LV), but much faster.

Usage

stepLVc(x0,t0,deltat,th=c(1,0.005,0.6))

Arguments

Value

A 2-vector representing the new state of the LV system.

See Also

StepGillespie, spnModels, simTs, simSample

Examples

# load the LV model
data(spnModels)
# create a stepping function
stepLV = StepGillespie(LV)
# step the function
print(stepLV(c(x1=50,x2=100),0,1))
# simulate a realisation of the process and plot it
out = simTs(c(x1=50,x2=100),0,100,0.1,stepLV)
plot(out)
# now use "stepLVc" instead...
out = simTs(c(x1=50,x2=100),0,100,0.1,stepLVc)
plot(out)