## ----setup, include=FALSE------------------------------------------------------------------------- knitr::opts_chunk$set(collapse=TRUE, cache=FALSE) ## ----set-options, echo=FALSE---------------------------------------------------------------------- options(width=1e2) ## ----eval=FALSE, include=FALSE-------------------------------------------------------------------- # # We can use this to type 'rmd' in the console, and the markdown html document is then knitted fast. # render.rmd("rmd") ## ------------------------------------------------------------------------------------------------- ## Load libraries library(ctsmTMB) library(ggplot2) ## plots ## Create model model <- newModel() model$addSystem(dx ~ theta * (a*u^2-cos(2*pi*t*u) - x) * dt + sigma_x*dw) model$addObs(y ~ x) model$setVariance(y ~ sigma_y^2) model$addInput(u) ## Set parameter values ## note: not strictly necessary to set lower/upper bounds model$setParameter( theta = c(initial = 2, lower = 0, upper = 100), a = c(initial = 1, lower = 1e-5, upper = 100), sigma_x = c(initial = 0.2, lower = 1e-5, upper = 5), ## fix sigma_y to 0.05 by not giving any upper/lower bounds sigma_y = c(initial = 5e-2) ) ## Set initial state mean and covariance ## note: diag(1) is not strictly needed model$setInitialState(list(1, 1e-1*diag(1))) ## ------------------------------------------------------------------------------------------------- ## set true parameters, and create data true.pars <- c(theta=4, a=10, sigma_x=1, sigma_y=0.05) dt.sim <- 1e-3 t.sim <- seq(0, 1, by=dt.sim) ## seed for input creation set.seed(20) u.sim <- cumsum(rnorm(length(t.sim),sd=0.1)) # u.sim <- tanh(sin(u.sim)) df.sim <- data.frame(t=t.sim, y=NA, u=u.sim) ## set rng seeds for C++ states and observations cpp.seeds <- c(20,20) ## perform simulation sim <- model$simulate(data=df.sim, pars=true.pars, n.sims=1, silent=T, cpp.seeds = cpp.seeds) ## extract simulations y.sim <- sim$observations$y$i0[,1] ## extract observations as every 10th simulation iobs <- seq(1, length(t.sim), by=10) t.obs <- t.sim[iobs] y.obs <- y.sim[iobs] u.obs <- u.sim[iobs] ## create data df.obs <- data.frame( t = t.obs, u = u.obs, y = y.obs ) ## ----echo=TRUE, fig.height=5,fig.width=9,out.width="100%", fig.align='center'--------------------- ggplot() + geom_line(aes(x=t.sim, y=y.sim, color="y_t")) + geom_point(aes(x=t.obs, y=y.obs, color="y_k")) + geom_line(aes(x=t.sim, y=u.sim, color="u_t")) + geom_point(aes(x=t.obs, y=u.obs, color="u_k")) + scale_color_discrete( labels=expression(y[t],y[k],u[t],u[k]), breaks = c("y_t", "y_k", "u_t", "u_k"), ) + ctsmTMB:::getggplot2theme() + theme(legend.text = ggplot2::element_text(size=15)) + labs(color="",x="Time",y="") ## ------------------------------------------------------------------------------------------------- ## set the first observation to NA df.temp <- df.obs df.temp$y[1] <- NA ## predict pred <- model$predict(data=df.temp, pars=true.pars, k.ahead=nrow(df.temp)-1) ## ------------------------------------------------------------------------------------------------- head(pred$states) head(pred$observations) ## ----message=FALSE-------------------------------------------------------------------------------- ## create vectors with different theta's pars <- lapply(c(1,2,4,8), function(x) c(x, 10, 1, 0.05)) ## predict for the various theta values pred1 = model$predict(df.obs, pars=pars[[1]]) pred2 = model$predict(df.obs, pars=pars[[2]]) pred3 = model$predict(df.obs, pars=pars[[3]]) pred4 = model$predict(df.obs, pars=pars[[4]]) ## ----echo=TRUE, fig.height=5,fig.width=9,out.width="100%", fig.align='center'--------------------- plot.df <- data.frame(t = pred$states[,"t.j"], t.obs = t.obs, y.obs = y.obs, x1=pred1$states[,"x"], x2=pred2$states[,"x"], x3=pred3$states[,"x"], x4=pred4$states[,"x"]) ggplot(data=plot.df) + geom_line(aes(x=t, y=x1, color="x1")) + geom_line(aes(x=t, y=x2, color="x2")) + geom_line(aes(x=t, y=x3, color="x3")) + geom_line(aes(x=t, y=x4, color="x4")) + geom_point(aes(x=t.obs, y=y.obs, color="y")) + scale_color_discrete( breaks = c("x1", "x2", "x3", "x4", "y"), labels=expression(theta*"=1",theta*"=10",theta*"=50",theta*"=100",y[k]) ) + ctsmTMB:::getggplot2theme() + theme(legend.text = ggplot2::element_text(size=15)) + labs(color="",x="Time",y="") ## ------------------------------------------------------------------------------------------------- fit = model$estimate(df.obs) print(fit) ## ------------------------------------------------------------------------------------------------- pred.horizon <- 15 pred = model$predict(df.obs, k.ahead=pred.horizon) ## ------------------------------------------------------------------------------------------------- pred.obs = pred$states[pred$states[,"k.ahead"]==pred.horizon,] ## ----echo=FALSE, fig.height=5,fig.width=9,out.width="100%", fig.align='center', fig.alt='The plot shows 25-step predictions with associated 95% prediction interval in grey, against the observations.'---- t.j <- pred.obs[,"t.j"] x <- pred.obs[,"x"] x.sd <- sqrt(pred.obs[,"var.x"]) ggplot() + geom_line(aes(x=t.j, y=x,color="25-Step Predictions")) + geom_ribbon(aes(x=t.j, ymin=x-2*x.sd, ymax=x + 2*x.sd),fill="grey",alpha=0.5) + geom_point(aes(x=t.obs, y.obs,color="Observations")) + labs(color="",x="Time",y="") + ctsmTMB:::getggplot2theme() ## ------------------------------------------------------------------------------------------------- pdf <- data.frame(pred$states, y.data = pred$observations[,"y.data"]) rmse <- data.frame( k.ahead = 0:pred.horizon, rmse = sapply( with(pdf, split(pdf, k.ahead)), function(df) { sqrt(mean((df$x - df$y.data)^2, na.rm = TRUE)) }) ) ## ----echo=FALSE, fig.height=5,fig.width=9,out.width="100%", fig.align='center', fig.alt='The plot shows the Root Mean Square Error score for prediction horizons from 0 to 25.'---- ggplot(data=rmse) + geom_line(aes(k.ahead, rmse), color="steelblue") + geom_point(aes(k.ahead, rmse), color="red") + labs( title = "Root-Mean Square Errors for Different Prediction Horizons", x = "Prediction Steps", y = "Root-Mean-Square Errors" ) + ctsmTMB:::getggplot2theme() ## ----eval=FALSE----------------------------------------------------------------------------------- # model$predict(data, # pars = NULL, # method = "ekf", # ode.solver = "rk4", # ode.timestep = diff(data$t), # k.ahead = nrow(data)-1, # return.k.ahead = 0:min(k.ahead, nrow(data)-1), # return.variance = c("marginal", "none", "covariance", "correlation"), # ukf.hyperpars = c(1, 0, 3), # initial.state = self$getInitialState(), # estimate.initial.state = private$estimate.initial, # use.cpp = TRUE, # silent = FALSE, # ...)