Nothing
## ----opts, echo=FALSE, results='hide'-----------------------------------------
knitr::opts_chunk$set(dev='png', dpi=150)
## ----install1, eval=FALSE-----------------------------------------------------
# install.packages("deBInfer")
## ----install2, eval=FALSE-----------------------------------------------------
# if (require("devtools")){
# #install deBInfer from github
# devtools::install_github("pboesu/debinfer")
# }
## ----loadlib, message=FALSE---------------------------------------------------
library(deBInfer)
## ----ode-def, message=FALSE,warning=FALSE-------------------------------------
library(deSolve)
logistic_model <- function (time, y, parms) {
with(as.list(c(y, parms)), {
dN <- r * N * (1 - N / K)
list(dN)
})
}
## ----pars-def-----------------------------------------------------------------
y <- c(N = 0.1)
parms <- c(r = 0.1, K = 10)
times <- seq(0, 120, 1)
out <- ode(y, times, logistic_model, parms, method="lsoda")
## ---- echo=FALSE--------------------------------------------------------------
plot(out)
## -----------------------------------------------------------------------------
set.seed(143)
#force include the first time-point (t=0)
N_obs <- as.data.frame(out[c(1,runif(35, 0, nrow(out))),])
## -----------------------------------------------------------------------------
# add lognormal noise
parms["sdlog.N"] <- 0.05
N_obs$N_noisy <- rlnorm(nrow(N_obs), log(N_obs$N),parms["sdlog.N"])
#observations must be ordered for solver to work
N_obs <- N_obs[order(N_obs$time),]
#Plot true and noisy observations
plot(N_obs$time, N_obs$N, ylim=c(0, max(N_obs$N,N_obs$N_noisy)))
points(N_obs$time, N_obs$N_noisy, col="red")
## ----obs-model----------------------------------------------------------------
# the observation model
logistic_obs_model <- function(data, sim.data, samp){
llik.N <- sum(dlnorm(data$N_noisy, meanlog = log(sim.data[,"N"] + 1e-6),
sdlog = samp[["sdlog.N"]], log = TRUE))
return(llik.N)
}
## ----pars, results="hide", message=FALSE--------------------------------------
library(deBInfer)
r <- debinfer_par(name = "r", var.type = "de", fixed = FALSE,
value = 0.5, prior = "norm", hypers = list(mean = 0, sd = 1),
prop.var = 0.0001, samp.type="rw")
K <- debinfer_par(name = "K", var.type = "de", fixed = FALSE,
value = 5, prior = "lnorm", hypers = list(meanlog = 1, sdlog = 1),
prop.var = 0.1, samp.type = "rw")
sdlog.N <- debinfer_par(name = "sdlog.N", var.type = "obs", fixed = FALSE,
value = 0.05, prior = "lnorm", hypers = list(meanlog = 0, sdlog = 1),
prop.var = c(3,4), samp.type = "rw-unif")
## ----inits--------------------------------------------------------------------
N <- debinfer_par(name = "N", var.type = "init", fixed = TRUE, value = 0.1)
## ----setup--------------------------------------------------------------------
mcmc.pars <- setup_debinfer(r, K, sdlog.N, N)
## ----deBinfer, results="hide", message=FALSE----------------------------------
# do inference with deBInfer
# MCMC iterations
iter <- 5000
# inference call
mcmc_samples <- de_mcmc(N = iter, data = N_obs, de.model = logistic_model,
obs.model = logistic_obs_model, all.params = mcmc.pars,
Tmax = max(N_obs$time), data.times = N_obs$time, cnt = 500,
plot = FALSE, verbose.mcmc = FALSE, solver = "ode")
## ----message=FALSE, warning=FALSE,fig.width = 8, fig.height = 8---------------
plot(mcmc_samples)
## -----------------------------------------------------------------------------
burnin <- 250
pairs(mcmc_samples, burnin = burnin, scatter=TRUE, trend=TRUE)
## ----post-dens, fig.height=5--------------------------------------------------
par(mfrow = c(1,3))
post_prior_densplot(mcmc_samples, burnin = burnin, param = "r")
abline(v = 0.1, col = "red", lty = 2)
post_prior_densplot(mcmc_samples, burnin = burnin, param = "K")
abline(v = 10, col = "red", lty = 2)
post_prior_densplot(mcmc_samples, burnin = burnin, param = "sdlog.N")
abline(v = 0.05, col = "red", lty = 2)
## ----post-sims----------------------------------------------------------------
post_traj <- post_sim(mcmc_samples, n=500, times=0:100, burnin=burnin, output = "all", prob = 0.95)
## ----post-sims-plot-----------------------------------------------------------
#median and HDI
plot(post_traj, plot.type = "medianHDI", lty = c(2,1), lwd = 3, col = c("red","grey20"),
panel.first = lines(out, col = "darkblue", lwd = 2)
)
legend("topleft", legend = c("posterior median", "95% HDI", "true model"),
lty = c(2,1,1), lwd = c(3,2,2), col = c("red","grey20","darkblue"),
bty = "n")
## ----post-sims-ensemble-------------------------------------------------------
plot(post_traj, plot.type = "ensemble", col = "#FF000040")
## ----epsilon-sims, eval=FALSE-------------------------------------------------
# #reformulate the observation model such that epsilon is an observation parameter
# logistic_obs_model_eps <- function(data, sim.data, samp){
#
# llik.N <- sum(dlnorm(data$N_noisy, meanlog = log(sim.data[,"N"] + samp[["epsilon"]]),
# sdlog = samp[["sdlog.N"]], log = TRUE))
# return(llik.N)
# }
#
# #declare epsilon
# epsilon <- debinfer_par(name = "epsilon", var.type = "obs", fixed = TRUE,
# value = 1e-6)
#
# #set up MCMC parameters
# mcmc.pars <- setup_debinfer(r, K, sdlog.N, epsilon, N)
## ----epsilon-sims-ctd, eval=FALSE---------------------------------------------
# # define a range of epsilon values
# eps_range <- c(1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1)
#
# # conduct inference for each value of eps_range
# eps_sims <- lapply(eps_range, function(x){
# #reassign value of epsilon
# mcmc.pars$epsilon$value <- x
# eps_samples <- de_mcmc(N = 10000, data = N_obs, de.model = logistic_model,
# obs.model = logistic_obs_model_eps, all.params = mcmc.pars,
# Tmax = max(N_obs$time), data.times = N_obs$time, cnt = 500,
# plot = FALSE, solver = "ode")
# return(eps_samples)
# })
#
# # add names to the inference runs
# names(eps_sims) <- as.character(eps_range)
#
# # collate the samples in a data.frame for easy plotting while omitting
# # a burnin of 1000 samples
# eps_samps <- plyr::ldply(eps_sims, function(x)(x$samples[1000:10000,]), .id = "epsilon")
#
# # plot sensitivity analysis
# par(mfrow = c(1,3), mar = c(3,3,2,2), mgp = c(1.7,0.5,0))
# beanplot::beanplot(eps_samps$r ~ eps_samps$epsilon, log = "", what = c(0,1,1,0),
# col = "darkgrey", bw = 2e-4, maxwidth = 0.6, xlab = "epsilon",
# ylab = "r")
# legend("bottomleft", legend = c("posterior mean", "true value"), lty = c(1,3),
# col = c("black", "red"), lwd = 2, bty = "n")
# abline(h = 0.1, col = "red", lwd = 2, lty = 3)
# beanplot::beanplot(eps_samps$K ~ eps_samps$epsilon, log = "", what = c(0,1,1,0),
# col = "darkgrey", bw = 2e-2, maxwidth = 0.6, xlab = "epsilon"
# , ylab = "K")
# abline(h = 10, col = "red", lwd = 2, lty = 3)
# beanplot::beanplot(eps_samps$sdlog.N ~ eps_samps$epsilon, log="", what=c(0,1,1,0),
# col = "darkgrey", bw = 2e-3, maxwidth = 0.6, xlab = "epsilon",
# ylab = "sdlog.N")
# abline(h = 0.05, col = "red", lwd = 2, lty = 3)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.