stan_ode: Simulating/Fitting ODE Systems via Stan

In imadmali/stanode: R wrapper for modeling ODEs in Stan

Description

Ordinary differential equation (ODE) simulation and parameter estimation. The behavior is similar to the ode function in the deSolve package. Particularly, the user may use the same ODE function, inital state values, parameter values used in ode.

Usage

stan_ode(func, state, pars, times, t0 = NULL, integrator = c("rk45", "bdf"),
  events = NULL, sampling = FALSE, y = NULL, likelihoods = NULL,
  priors = NULL, ...)

Arguments

func	A function (currently if-else statements are not supported).
state	A named vector of the initial conditions of the state variables.
pars	A named vector of the parameter values.
times	A sequence of time steps.
t0	Initial time. The default (NULL) is t0 = times[1] - 1e-6
integrator	The type of integrator to use.
events	A dataframe that defines the event schedule in the ODE system. See details below.
sampling	A logical declaration as to whether you want sample the parameters or initial state values. The default is FALSE.
...	Optional parameters for stan.

Details

Currently the user cannot use the event dataframe used in ode. The dimensions of the events dataframe equals [number of events] by [number of state variables + 2]. There must be a column denoting the time of each event, the method to apply to the state variables, and the value of each state variable that the user wants to apply to the last simulated value using the aforementioned method.

Some examples of event schedules are provided below. The event schedule below is adding 5 to the last simulated state variable y1 in all events accept for the one at taking place time = 40 where 5 is being added to y1 and 2 is being added to y2.

time	y1	y2	...	method
10	5	0	...	add
20	5	0	...	add
30	5	0	...	add
40	5	2	...	add

The event schedule below is multiplying 1.5 with the last simulated value of y1 and leaving y2 unchanged.

time	y1	y2	...	method
10	1.5	1	...	multiply
20	1.5	1	...	multiply
30	1.5	1	...	multiply

Value

A list that contains the simulations and the stanfit object.

See Also

ode, stan.

Examples

## Not run: 
# Simple ODE
f <- function(y, t, p) {
  dy1 <- y2
  dy2 <- -y1 - theta1 * y2
  return(dy1 = dy1, dy2 = dy2)
}
fit <- stan_ode(f, state = c("y1" = 2, "y2" = 5),
                pars = c("theta1" = 0.5),
                times = seq(1,10,by=0.01), t0 = 0,
                integrator = "bdf",
                sampling = FALSE)
sims <- extract(fit, pars = "y_hat")
sims <- unname(unlist(sims))

# Simple Harmonic Oscillator
sho <- function(y, p ,t) {
  dy1 = y2
  dy2 = -y1 - theta * y2
  return(list(dy1, dy2))
}

fit <- stan_ode(sho, state = c("y1" = 1, "y2" = 0),
                pars = c("theta" = 0.15),
                times = seq(1,50,by=0.1), t0 = 0,
                integrator = "bdf",
                sampling = FALSE)

sims <- extract(fit, pars = "y_hat")$y_hat[1,,]
plot(sims[,1], sims[,2], type = "l", lwd = 2,
     xlab = "y1", ylab = "y2", main = "Simple Harmonic Oscillator")

## End(Not run)

imadmali/stanode documentation built on May 3, 2019, 11:48 p.m.