StepSDE: Create a function for advancing the state of an SDE model by...

In smfsb: Stochastic Modelling for Systems Biology

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

drift	A function representing the drift vector of the SDE model (corresponding roughly to the RHS of an ODE model). drift should have prototype drift(x,t,...), with x representing current system state and t representing current system time. The value of the function should be a vector of the same dimension as x, representing the infinitesimal mean of the Ito SDE.
diffusion	A function representing the diffusion matrix of the SDE model (the square root of the infinitesimal variance matrix). diffusion should have prototype diffusion(x,t,...), with x representing current system state and t representing current system time. The value of the function should be a square matrix with both dimensions the same as the length of x.
dt	Time step to be used by the simple Euler-Maruyama integration method. Defaults to 0.01.

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)

smfsb documentation built on May 2, 2019, 5:13 a.m.