require(deSolve) require(rstan) ### simulate data m <- function(t, y, p) { with(as.list(c(y, p)), { dydt <- y * r * (1 - y/k) return(list(c(dydt))) }) } sim_ode <- as.data.frame( rk(func = m, y = c(y=1), parms = c(r=0.2, k = 200), times = seq(0,100,1), method="rk45dp6") ) # sim_ode <- as.data.frame( # ode(func = m, y = c(y=1), parms = c(r=0.2, k = 50), times = seq(0,100,1), method="iteration") # ) plot(sim_ode$y, type="l") y_noise <- rlnorm(nrow(sim_ode), log(sim_ode$y), 0.1) stan_y <- y_noise[-1] # estimate parameters using the Stan ode model fit_stan_ode <- stan( file = "logistic-ode.stan", init = 0, data = list( y = stan_y, T = length(stan_y), ts = seq(1, length(stan_y), 1) ) ) print(fit_stan_ode, pars = c("y_init","sigma","theta")) ps_ode <- as.matrix(fit_stan_ode) plot(sim_ode$y, type="n", ylim=c(0,250), main = "ODE") for(i in sample(1:nrow(ps_ode), 30, replace = TRUE)){ lines(c(ps_ode[i,grep("y_init", colnames(ps_ode))], ps_ode[i,grep("\\bz\\b", colnames(ps_ode))]), lwd = 0.5) } points(y_noise, pch = 16, cex = 1, col = "red") lines(sim_ode$y, lwd=3, col="blue") # estimate parameters using the Stan dts model fit_stan_dts <- stan( file = "logistic-dts.stan", init = 0, data = list( y = stan_y, T = length(stan_y) ) ) print(fit_stan_dts, pars = c("y_init","sigma","theta")) ps_dts <- as.matrix(fit_stan_dts) plot(sim_ode$y, type="n", ylim=c(0,250), main = "DTS") for(i in sample(1:nrow(ps_dts), 30, replace = TRUE)){ lines(c(ps_dts[i,grep("y_init", colnames(ps_dts))], ps_dts[i,grep("\\bz\\b", colnames(ps_dts))]), lwd = 0.5) } points(y_noise, pch = 16, cex = 1, col = "red") lines(sim_ode$y, lwd=3, col="blue")