context("Testing inference procedure using a simple logistic model with joint proposals")
test_that("Inference on simulated data with known inits. ", {
#testthat::skip_on_cran()
# define logistic model
logistic_model <- function(time, y, parms) {
with(as.list(c(y, parms)), {
dN <- r * N * (1 - N / K)
list(dN)
})
}
testthat::expect_true(exists("logistic_model"))
# set initial value for simulation
y <- c(N = 0.1)
# set parameter values
parms <- c(r = 0.1, K = 10)
# set simulation time points
times <- seq(0, 120, 1)
# solve ODE
out <- ode(y, times, logistic_model, parms, method='lsoda')
# sample from simulated data
set.seed(143)
N_obs <- as.data.frame(out[c(1,runif(35, 0, nrow(out))),]) #force include the first time-point (t=0)
# add lognormal noise
parms['logsd.N'] <- 0.01
N_obs$N_noisy <- rlnorm(nrow(N_obs), log(N_obs$N),(parms['logsd.N']))
# observations must be ordered for solver to work
N_obs <- N_obs[order(N_obs$time),]
# define an observation model
# NB: lognormal errors are not great really for the obs model - should be changed to sth that actually allows N to be zero instead of using epsilon correction
logistic_obs_model<-function(data, sim.data, samp){
llik.N<-sum(dlnorm(data$N_noisy, meanlog=log(sim.data[,"N"]+1e-6), sdlog=samp[['logsd.N']], log=TRUE))
llik<-llik.N
return(llik)
}
r <- debinfer_par(name = "r", var.type = "de", fixed = FALSE,
value = 0.5, prior="norm", hypers=list(mean = 0, sd = 1),
prop.var=5e-5, samp.type="rw", joint="sigma1")
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", joint="sigma1")
#define covariance matrix
sigma1 <- debinfer_cov(c("r", "K"), sigma = matrix(c(5e-05, -1e-05,-1e-05, 1e-03),2, byrow = TRUE), name = "sigma1")
logsd.N <- debinfer_par(name = "logsd.N", var.type = "obs", fixed = FALSE,
value = 1, prior="lnorm", hypers=list(meanlog = 0, sdlog = 1),
prop.var=c(1,2), samp.type="rw-unif")
#we also need to provide an initial condition for the DE
N <- debinfer_par(name = "N", var.type = "init", fixed = TRUE,
value = 0.1)
mcmc.pars <- setup_debinfer(r, K, logsd.N, N, sigma1)
# do inference with deBInfer
# MCMC iterations
iter = 7000
# define burnin
burnin = 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=iter+1,
plot=FALSE, sizestep=0.1, solver=1)
#add more tests here checking the integrity & contents of the returned data structure
#check accuracy of estimation (threshold is 5% of true de parameter value and 20% of true observation noise)
expect_equal(unname(mean(mcmc_samples$samples[burnin:iter,"r"])/parms["r"]),1,tolerance = 5e-2)
expect_equal(unname(mean(mcmc_samples$samples[burnin:iter,"K"])/parms["K"]),1,tolerance = 5e-2)
expect_equal(unname(mean(mcmc_samples$samples[burnin:iter,"logsd.N"])/parms["logsd.N"]),1,tolerance = 2e-1)
#test utility function for checking results class
expect_equal(is.debinfer_result(mcmc_samples), TRUE)
#test extractor functions
expect_equal(deinits(mcmc_samples), c(N=0.1))
expect_equal(depars(mcmc_samples), c(r=0.5, K=5))
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.