diff --git a/inst/doc/chytrid_dede_example.R b/inst/doc/chytrid_dede_example.R
index b7d1597..ae57b33 100644
--- a/inst/doc/chytrid_dede_example.R
+++ b/inst/doc/chytrid_dede_example.R
@@ -113,7 +113,7 @@ iter = 500
 # inference call
 dede_rev <- de_mcmc(N = iter, data=chytrid, de.model=CSZ.dede,
                                obs.model=chytrid_obs_model, all.params=mcmc.pars,
-                               Tmax = max(chytrid$time), data.times=c(0,chytrid$time), cnt=iter %/% 10,
+                               Tmax = max(chytrid$time), data.times=c(0,chytrid$time), cnt=50,
                                plot=FALSE, sizestep=0.1, solver="dede", verbose = TRUE)
 
 ## ----message=FALSE, warning=FALSE,fig.width = 8, fig.height = 8----------
@@ -121,19 +121,43 @@ par(mfrow = c(3,4))
 plot(dede_rev, ask=FALSE, auto.layout=FALSE)
 
 ## ---- fig.width = 8, fig.height = 8--------------------------------------
-burnin = 50
+burnin = 100
 pairs(dede_rev, burnin = burnin, scatter=TRUE, trend=TRUE)
 post_prior_densplot(dede_rev, burnin = burnin)
-summary(dede_rev)
+
+## ------------------------------------------------------------------------
+par(mfrow=c(2,3), mgp=c(2.2, 0.8, 0))
+#define a fancy y axis label
+ylabel = expression(paste(Pr,"(", theta,"|", "Y", ")"))
+#plot the individual parameters
+post_prior_densplot(dede_rev, param="sr",xlab=expression(theta), 
+                    ylab=ylabel, show.obs=FALSE, xlim=c(0,8), 
+                    main=expression(paste("s",phantom()[{paste("r")}])))
+legend("topright", legend=c("Posterior","Prior"), lty = 1, col = c("black", "red"))
+post_prior_densplot(dede_rev, param="fs",xlab=expression(theta), 
+                    ylab=ylabel, show.obs=FALSE, xlim=c(-0.1,1.1), 
+                    main=expression(paste("f",phantom()[{paste("s")}])))
+post_prior_densplot(dede_rev, param="ds",xlab=expression(theta),
+                    ylab=ylabel, show.obs=FALSE, xlim=c(0,3), 
+                    main=expression(paste("d",phantom()[{paste("s")}])))
+post_prior_densplot(dede_rev, param="muz",xlab=expression(theta),
+                    ylab=ylabel, show.obs=FALSE, xlim=c(0,6), 
+                    main=expression(paste(mu,phantom()[{paste("Z")}])))
+post_prior_densplot(dede_rev, param="eta",xlab=expression(theta),
+                    ylab=ylabel, show.obs=FALSE, xlim=c(0,50), ylim=c(0,0.2), 
+                    main=expression(eta))
+post_prior_densplot(dede_rev, param="Tmin",xlab=expression(theta),
+                    ylab=ylabel, show.obs=FALSE, xlim=c(1.5,6.5), 
+                    main=expression(paste("T",phantom()[{paste("min")}])))
 
 ## ----post-sims-----------------------------------------------------------
-post_traj <- post_sim(dede_rev, n=20, times=0:100, burnin=burnin, output = 'all')
+post_traj <- post_sim(dede_rev, n=100, times=0:100, burnin=burnin, output = 'all', prob = 0.95)
 
 ## ----post-sims-plot, fig.width = 8, fig.height = 3-----------------------
 #median and HDI
 par(mfrow=c(1,3))
 plot(post_traj, plot.type = "medianHDI", auto.layout = FALSE)
-legend("topright", legend=c("posterior median", "95% HDI"), lty=1, col=c("red","grey"))
+legend("topright", legend=c("posterior median", "95% HDI"), lty=1, col=c("red","grey"), bty='n')
 
 
 ## ----post-sims-ensemble, fig.width = 8, fig.height = 6-------------------
diff --git a/inst/doc/chytrid_dede_example.Rmd b/inst/doc/chytrid_dede_example.Rmd
index f5e93e5..1f78341 100644
--- a/inst/doc/chytrid_dede_example.Rmd
+++ b/inst/doc/chytrid_dede_example.Rmd
@@ -4,7 +4,6 @@ author: "Philipp H Boersch-Supan and Leah R Johnson"
 date: '`r format(Sys.Date(), "%B %d, %Y")`'
 output: 
   pdf_document: 
-    keep_tex: yes
     number_sections: yes
 bibliography: debinfer.bib
 vignette: >
@@ -15,7 +14,9 @@ vignette: >
 
 #Preliminaries
 
-This examples assumes that `deBInfer` is installed and loaded. If this is not the case it needs to be installed from github, which requires the devtools package.
+This vignette illustrates the steps needed to perform inference for a DDE model of population growth in a fungal pathogen. A detailed description of the rationale behind the Bayesian inference approach for differential equations can be found in the preprint describing the `deBInfer` package [@deBInfer]. 
+
+This example assumes that `deBInfer` is installed and loaded. If this is not the case it needs to be installed from github, which requires the devtools package [@devtools].
 ```{r install1, eval=FALSE}
 install.packages("devtools")
 ```
@@ -107,9 +108,7 @@ The log-likelihood of the data given the parameters, underlying model, and initi
 \ell(\mathbf{Z}|\mathbf{\theta})=\sum\limits^n_{t}Z_t \text{ log }\lambda-n\lambda
 \end{equation}
 
-**Explain epsilon correction**
-
-This can be translated into an observation model function for `deBInfer`. The observation model function must have three named arguments `data`, `sim.data`, and `samp`, as these are used by the MCMC procedure to pass in the data (as a `data.frame`), the current state of the Markov chain (as a named vector), and the associated DE model solution (as a matrix-like object of class `deSolve`). We can access these inputs to define the data likelihood. In this case we have repeat measurements for each time point, so we iterate over the unique timepoints in `data$time`, and then calculate the sum log-likelihood over all matching `data$count` observations using the current value of the state variable $Z$ from the DE model solution at this point in the Markov chain.
+This can be translated into an observation model function for `deBInfer`. The observation model function must have three named arguments `data`, `sim.data`, and `samp`, as these are used by the MCMC procedure to pass in the data (as a `data.frame`, i.e. indexed using the `$` operator and the column names of the input data), the current state of the Markov chain (as a named vector), and the associated DE model solution (as a matrix-like object of class `deSolve`, i.e. indexed using the `[ ]` operator and the declared names of the state variables). We can access these inputs to define the data likelihood. In this case we have repeat measurements for each time point, so we iterate over the unique timepoints in `data$time`, and then calculate the sum log-likelihood over all matching `data$count` observations using the current value of the state variable $Z$ from the DE model solution at this point in the Markov chain.
 
 ```{r obs-model}
 # observation model
@@ -126,8 +125,11 @@ chytrid_obs_model<-function(data, sim.data, samp){
   return(llik)
 }
 ```
+Following Johnson et al. [-@johnson2013bayesian] we employ a small correction `ec` that is needed because the DE solution can equal zero, whereas the parameter `lambda` of the Poisson likelihood must be strictly positive.
 
 
+## Parameter declaration
+
 We continue by defining the parameters for inference
 ```{r vars}
 sr <- debinfer_par(name = "sr", var.type = "de", fixed = FALSE,
@@ -161,14 +163,16 @@ S <- debinfer_par(name = "S", var.type = "init", fixed = TRUE, value = 0)
 Z <- debinfer_par(name = "Z", var.type = "init", fixed = TRUE, value = 0)
 ```
 
-**Explain importance of init order**
+## MCMC Inference 
+
+The declared parameters are then collated using the `setup_debinfer` function. Note that **the initial values must be entered in the same order, as they are specified in the DE model function**, as the solver matches these values by position, rather than by name. More details can be found in `?deSolve::dede`. The remaining parameters can be entered in any order.
 
 ```{r mcmc-setup}
 # ----setup---------------------------------------------------------------
 mcmc.pars <- setup_debinfer(sr, fs, ds, muz, eta, Tmin, C, S, Z)
 ```
 
-And then do the actual inference
+`de_mcmc` is the workhorse of the package and runs the MCMC estimation. The progress of the MCMC procedure can be monitored using the `cnt`, `plot` and `verbose` options: Every `cnt` iterations the function will print out information about the current state, and, if `plot=TRUE`, traceplots of the chains will be plotted. Setting `verbose=TRUE` will print additional information. Note that frequent plotting will substantially slow down the MCMC sampler, and should be used only on short runs when tuning the sampler.
 ```{r de_mcmc, results='hide'}
 # do inference with deBInfer
 # MCMC iterations
@@ -176,37 +180,67 @@ iter = 500
 # inference call
 dede_rev <- de_mcmc(N = iter, data=chytrid, de.model=CSZ.dede,
                                obs.model=chytrid_obs_model, all.params=mcmc.pars,
-                               Tmax = max(chytrid$time), data.times=c(0,chytrid$time), cnt=iter %/% 10,
+                               Tmax = max(chytrid$time), data.times=c(0,chytrid$time), cnt=50,
                                plot=FALSE, sizestep=0.1, solver="dede", verbose = TRUE)
 ```
+Note that the number of iterations was set to 500 to keep the build time of the vignette within acceptable limits. For the inference results reported in Boersch-Supan and Johnson [-@deBInfer] the MCMC procedure was run for 100000 iterations, which took about 4 hours on a 2014 Apple Mac mini with a 2.6 GHz Intel i5 processor.
+
+## MCMC diagnostics
 
 We plot and summarize the MCMC chains
 ```{r message=FALSE, warning=FALSE,fig.width = 8, fig.height = 8}
 par(mfrow = c(3,4))
 plot(dede_rev, ask=FALSE, auto.layout=FALSE)
 ```
-From the traceplot we can see that the burnin period is **about 200 samples**. We can remove the burnin and have a look at parameter correlations, and the overlap between the posterior and prior densities.
+From the traceplot we can see that with only 500 iterations the chains have neither mixed well, nor reached stationarity. For demonstration purposes we remove a burnin period of 100 samples and have a look at parameter correlations, and the overlap between the posterior and prior densities.
 ```{r, fig.width = 8, fig.height = 8}
-burnin = 50
+burnin = 100
 pairs(dede_rev, burnin = burnin, scatter=TRUE, trend=TRUE)
 post_prior_densplot(dede_rev, burnin = burnin)
-summary(dede_rev)
 ```
-We simulate DE model trajectories from the posterior and calculate the HPD interval for the deterministic part of the model.
+More control over the `post_prior_densplot` can be achieved by plotting the parameters individually, using the `param` option. This way the x and y limits of the plots can be adjusted to show a larger portion of the prior support, and fancy labels can be added.
+```{r}
+par(mfrow=c(2,3), mgp=c(2.2, 0.8, 0))
+#define a fancy y axis label
+ylabel = expression(paste(Pr,"(", theta,"|", "Y", ")"))
+#plot the individual parameters
+post_prior_densplot(dede_rev, param="sr",xlab=expression(theta), 
+                    ylab=ylabel, show.obs=FALSE, xlim=c(0,8), 
+                    main=expression(paste("s",phantom()[{paste("r")}])))
+legend("topright", legend=c("Posterior","Prior"), lty = 1, col = c("black", "red"))
+post_prior_densplot(dede_rev, param="fs",xlab=expression(theta), 
+                    ylab=ylabel, show.obs=FALSE, xlim=c(-0.1,1.1), 
+                    main=expression(paste("f",phantom()[{paste("s")}])))
+post_prior_densplot(dede_rev, param="ds",xlab=expression(theta),
+                    ylab=ylabel, show.obs=FALSE, xlim=c(0,3), 
+                    main=expression(paste("d",phantom()[{paste("s")}])))
+post_prior_densplot(dede_rev, param="muz",xlab=expression(theta),
+                    ylab=ylabel, show.obs=FALSE, xlim=c(0,6), 
+                    main=expression(paste(mu,phantom()[{paste("Z")}])))
+post_prior_densplot(dede_rev, param="eta",xlab=expression(theta),
+                    ylab=ylabel, show.obs=FALSE, xlim=c(0,50), ylim=c(0,0.2), 
+                    main=expression(eta))
+post_prior_densplot(dede_rev, param="Tmin",xlab=expression(theta),
+                    ylab=ylabel, show.obs=FALSE, xlim=c(1.5,6.5), 
+                    main=expression(paste("T",phantom()[{paste("min")}])))
+```
+
+## Simulating posterior trajectories
+We simulate 100 DE model trajectories from the posterior and calculate the 95% highest posterior density interval for the deterministic part of the model.
 
 ```{r post-sims}
-post_traj <- post_sim(dede_rev, n=20, times=0:100, burnin=burnin, output = 'all')
+post_traj <- post_sim(dede_rev, n=100, times=0:100, burnin=burnin, output = 'all', prob = 0.95)
 ```
 
-We can visualise the median posterior trajectory and the associated highest posterior density interval using `r plot.type="medianHDI"` 
+We can visualise the median posterior trajectory and the highest posterior density interval using `r plot.type="medianHDI"` 
 ```{r post-sims-plot, fig.width = 8, fig.height = 3}
 #median and HDI
 par(mfrow=c(1,3))
 plot(post_traj, plot.type = "medianHDI", auto.layout = FALSE)
-legend("topright", legend=c("posterior median", "95% HDI"), lty=1, col=c("red","grey"))
+legend("topright", legend=c("posterior median", "95% HDI"), lty=1, col=c("red","grey"), bty='n')
 
 ```
-Alternatively we can plot an ensemble of posterior trajectories
+Alternatively we can plot an ensemble of posterior trajectories using `r plot.type="ensemble"` 
 ```{r post-sims-ensemble, fig.width = 8, fig.height = 6}
 plot(post_traj, plot.type = "ensemble", col = "#FF000040")
 ```
diff --git a/inst/doc/chytrid_dede_example.pdf b/inst/doc/chytrid_dede_example.pdf
index 7a0f1e6..f5d6c5e 100644
Binary files a/inst/doc/chytrid_dede_example.pdf and b/inst/doc/chytrid_dede_example.pdf differ
diff --git a/inst/doc/logistic_ode_example.R b/inst/doc/logistic_ode_example.R
index 6dee4bf..a2b125e 100644
--- a/inst/doc/logistic_ode_example.R
+++ b/inst/doc/logistic_ode_example.R
@@ -14,11 +14,13 @@ 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)
-})
+    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)
@@ -29,7 +31,8 @@ plot(out)
 
 ## ------------------------------------------------------------------------
 set.seed(143)
-N_obs <- as.data.frame(out[c(1,runif(35, 0, nrow(out))),]) #force include the first time-point (t=0)
+#force include the first time-point (t=0)
+N_obs <- as.data.frame(out[c(1,runif(35, 0, nrow(out))),]) 
 
 
 ## ------------------------------------------------------------------------
@@ -38,25 +41,17 @@ N_obs <- as.data.frame(out[c(1,runif(35, 0, nrow(out))),]) #force include the fi
   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")
 
-out_obs <- ode(y, c(0,N_obs$time), logistic_model, parms, method='lsoda')
-plot(out_obs, type="p")
-lines(out)
-
 ## ----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)
-                )
-
-    llik<-llik.N
-
-    return(llik)
+                       sdlog=samp[['sdlog.N']], log=TRUE))
+    return(llik.N)
   }
 
 
@@ -64,14 +59,14 @@ lines(out)
 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.005, samp.type="rw")
+                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.1, prior="lnorm", hypers=list(meanlog = 0, sdlog = 1),
+                value = 0.05, prior="lnorm", hypers=list(meanlog = 0, sdlog = 1),
                 prop.var=c(3,4), samp.type="rw-unif")
 
 ## ----inits---------------------------------------------------------------
@@ -87,26 +82,35 @@ mcmc.pars <- setup_debinfer(r, K, sdlog.N, N)
   # 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="ode")
+                          Tmax = max(N_obs$time), data.times=N_obs$time, cnt=500, 
+                          plot=FALSE, solver="ode")
 
 
 ## ----message=FALSE, warning=FALSE,fig.width = 8, fig.height = 8----------
 plot(mcmc_samples)
 
 ## ------------------------------------------------------------------------
-burnin = 200
+burnin = 250
 pairs(mcmc_samples, burnin = burnin, scatter=TRUE, trend=TRUE)
-post_prior_densplot(mcmc_samples, burnin = burnin)
-summary(mcmc_samples)
+
+## ------------------------------------------------------------------------
+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=200, times=0:100, burnin=burnin, output = 'all')
+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")
-legend("topleft", legend=c("posterior median", "95% HDI"), lty=1, col=c("red","grey"))
+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--------------------------------------------------
diff --git a/inst/doc/logistic_ode_example.Rmd b/inst/doc/logistic_ode_example.Rmd
index ff96a93..0b6eb7b 100644
--- a/inst/doc/logistic_ode_example.Rmd
+++ b/inst/doc/logistic_ode_example.Rmd
@@ -1,10 +1,9 @@
 ---
 title: "Bayesian inference for the logistic growth model"
-author: "Philipp H Boersch-Supan"
+author: "Philipp H Boersch-Supan and Leah R Johnson"
 date: '`r format(Sys.Date(), "%B %d, %Y")`'
 output: 
   pdf_document: 
-    keep_tex: yes
     number_sections: yes
 bibliography: debinfer.bib
 vignette: >
@@ -14,8 +13,9 @@ vignette: >
 ---
 
 #Preliminaries
+This vignette illustrates the steps needed to perform inference for an ODE model, the logistic growth model. A detailed description of the rationale behind the Bayesian inference approach for differential equations can be found in the preprint describing the `deBInfer` package [@deBInfer]. The vignette also illustrates how to simulate noisy observations from an existing DE model. Conducting inference on a model with known parameters is good practice, when developing models for actual data, as it provides validation of the inference procedure.
 
-This examples assumes that `deBInfer` is installed and loaded. If this is not the case it needs to be installed from github, which requires the devtools package.
+This examples assumes that `deBInfer` is installed and loaded. If this is not the case it needs to be installed from github, which requires the devtools package [@devtools].
 ```{r install1, eval=FALSE}
 install.packages("devtools")
 ```
@@ -34,97 +34,111 @@ library(deBInfer)
 
 #Defining the DE model
 
-We define a logistic growth ODE for deSolve.
+We define a logistic growth model 
+\begin{equation}\frac{dN}{dt}=rN\left(1-\frac{N}{K}\right).\end{equation}
+and implement it for `deSolve::ode` [@soetaert2010solving]:
 
 ```{r 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)
-})
+    with(as.list(c(y, parms)), {
+      dN <- r * N * (1 - N / K)
+      list(dN)
+      })
 }
+```
+
+We then pick a set of parameter values $(r=0.1, K=10)$ and an initial condition ($N=0.1$) and solve the model for these values.
+```{r 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')
 ```
 
-Which gives us the numerical solution
+We can plot this numerical solution.
 
 ```{r, echo=FALSE}
 plot(out)
 ```
 
-#Simulating observations
-Now we simulate a noisy dataset from this equation. We sample a random subset from the integration output 
+# Simulating observations
+Next we simulate a noisy dataset from this equation. We do this by sampling a random subset from the integration output 
 
 ```{r}
 set.seed(143)
-N_obs <- as.data.frame(out[c(1,runif(35, 0, nrow(out))),]) #force include the first time-point (t=0)
+#force include the first time-point (t=0)
+N_obs <- as.data.frame(out[c(1,runif(35, 0, nrow(out))),]) 
 
 ```
-and we "add" lognormal noise 
+
+and we "add" lognormal noise with a known "observation" standard deviation $\sigma^2_{obs}=0.05$.
 ```{r}
 # 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")
-
-out_obs <- ode(y, c(0,N_obs$time), logistic_model, parms, method='lsoda')
-plot(out_obs, type="p")
-lines(out)
 ```
 
 #Defining an observation model and parameters for inference
-We define an observation model. Note that we are sampling the log of the observation standard deviation, to ensure sampled values are strictly positive.
-We also use an epsilon correction for the meanlog, as the DE model can return values of 0 (or even less due to numerical precision). 
+
+The appropriate log-likelihood for these data takes the form 
+
+\begin{equation}
+\ell(\mathcal{Y}|\boldsymbol\theta) = \sum_t \ln\left(\frac{1}{\tilde{N_t}\sigma_{obs}\sqrt{2\pi}}\exp\left(-\frac{(\ln \tilde{N_t}-\ln (N_t+\varepsilon))^2}{2\sigma^2_{obs}}\right)\right)
+\end{equation}
+where $\tilde{N}_t$ are the observations, and $N_t$ are the predictions of the DE model given the current MCMC sample of the parameters $\boldsymbol{\theta}$. Further, $\varepsilon$  is a small correction needed, because the DE solution can equal zero (or less, depending on numerical precision) [@johnson2013bayesian]. We chose $\varepsilon= 10^{-6}$, which is on the same order as the  numerical precision of the default solver (`deSolve::ode`} with {`method = "lsoda"`).
+
+This can be translated into an observation model function for `deBInfer`. The observation model function must have three named arguments `data`, `sim.data`, and `samp`, as these are used by the MCMC procedure to pass in the data (as a `data.frame`, i.e. indexed using the `$` operator and the column names of the input data), the current state of the Markov chain (as a named vector), and the associated DE model solution (as a matrix-like object of class `deSolve`, i.e. indexed using the `[ ]` operator and the declared names of the state variables). We can access these inputs to define the data likelihood. 
+
+The user specifies the observation model such that it returns the summed log-likelihoods of the data. In this example the observations are in the data.frame column `N_noisy`, and the corresponding predicted states are in the column `N` of the matrix-like object`sim.data`
 ```{r 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)
-                )
-
-    llik<-llik.N
-
-    return(llik)
+                       sdlog=samp[['sdlog.N']], log=TRUE))
+    return(llik.N)
   }
 
 ```
 
-We declare the parameters for inference:
+All parameters that are used in the DE model and the observation model need to be declared for the inference procedure using the `debinfer_par()` function. The declaration describes the variable name, whether it is a DE or observation parameter and whether or not it is to be estimated. If the parameter is to be estimated, the user also needs to specify a prior distribution and a number of additional parameters for the MCMC procedure. `debinfer` currently supports priors from all probability distributions implemented in base R, as well as their truncated variants, as implemented in the `truncdist` package [@truncdist].
+
+We declare the DE model parameter $r$, assign a prior $r \sim \mathcal{N}(0,1)$ and a random walk sampler with a Normal kernel (`samp.type="rw"`) and  proposal variance of 0.005.
+Similarly, we declare $K\sim \ln\mathcal{N}(1,1)$ and $\ln(\sigma^2_{obs})\sim\mathcal{N}(0,1)$. Note that we are using the asymmetric uniform proposal distribution $\mathcal{U}(\frac{a}{b}\theta^{(k)}, \frac{b}{a}\theta^{(k)})$ for the variance parameter (`samp.type="rw-unif"`), as this ensures strictly positive proposals.
+
+
 ```{r 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.005, samp.type="rw")
+                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.1, prior="lnorm", hypers=list(meanlog = 0, sdlog = 1),
+                value = 0.05, prior="lnorm", hypers=list(meanlog = 0, sdlog = 1),
                 prop.var=c(3,4), samp.type="rw-unif")
 ```
 
-and we also need to provide an initial condition for the differential equation:
+Lastly, we provide an initial value $N_0=0.1$ for the DE:
 ```{r inits}
 N <- debinfer_par(name = "N", var.type = "init", fixed = TRUE, value = 0.1)
 ```
-All declared parameters are collated using the setup_debinfer function
+The declared parameters are then collated using the `setup_debinfer` function. Note that for models with more than one state variable, **the initial values must be entered in the same order, as they are specified in the DE model function**, as the solver matches these values by position, rather than by name. More details can be found in `?deSolve::ode`. The remaining parameters can be entered in any order.
 ```{r setup}
 mcmc.pars <- setup_debinfer(r, K, sdlog.N, N)
 ```
 
 #Conduct inference
-Finally we use deBInfer to estimate the parameters of the original model.
+Finally we use deBInfer to estimate the parameters of the original model. `de_mcmc` is the workhorse of the package and runs the MCMC estimation. The progress of the MCMC procedure can be monitored using the `cnt`, `plot` and `verbose` options: Every `cnt` iterations the function will print out information about the current state, and, if `plot=TRUE`, traceplots of the chains will be plotted. Setting `verbose=TRUE` will print additional information. Note that frequent plotting will substantially slow down the MCMC sampler, and should be used only on short runs when tuning the sampler.
 ```{r deBinfer, results="hide"}
 # do inference with deBInfer
   # MCMC iterations
@@ -132,8 +146,8 @@ Finally we use deBInfer to estimate the parameters of the original model.
   # 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="ode")
+                          Tmax = max(N_obs$time), data.times=N_obs$time, cnt=500, 
+                          plot=FALSE, solver="ode")
 
 ```
 
@@ -141,27 +155,47 @@ We plot and summarize the MCMC chains
 ```{r message=FALSE, warning=FALSE,fig.width = 8, fig.height = 8}
 plot(mcmc_samples)
 ```
-From the traceplot we can see that the burnin period is about 200 samples. We can remove the burnin and have a look at parameter correlations, and the overlap between the posterior and prior densities.
+
+From the traceplot we can see that the burnin period is about 250 samples. We can remove the burnin and have a look at parameter correlations using `pairs.debinfer_result`.
 ```{r}
-burnin = 200
+burnin = 250
 pairs(mcmc_samples, burnin = burnin, scatter=TRUE, trend=TRUE)
-post_prior_densplot(mcmc_samples, burnin = burnin)
-summary(mcmc_samples)
 ```
-We simulate DE model trajectories from the posterior and calculate the HPD interval for the deterministic part of the model.
+
+And we can look at the overlap between the posterior and prior densities using `post_prior_densplot`. 
+```{r}
+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)
+```
+Overlay the true underlying parameter values as vertical dashed lines shows that our estimation is recovering the true parameters.
+
+# Posterior simulations
+We then simulate DE model trajectories from the posterior and calculate the 95% highest posterior density interval for the deterministic part of the model.
 
 ```{r post-sims}
-post_traj <- post_sim(mcmc_samples, n=200, times=0:100, burnin=burnin, output = 'all')
+post_traj <- post_sim(mcmc_samples, n=500, times=0:100, burnin=burnin, output = 'all', prob = 0.95)
 ```
 
-We can visualise the median posterior trajectory and the associated highest posterior density interval using `r plot.type="medianHDI"` 
+We can visualise the median posterior trajectory and the associated highest posterior density interval using `plot.type="medianHDI"`. The `panel.first` option can be used to add the trajectory of the true model, which is stored in `out`.
 ```{r post-sims-plot}
 #median and HDI
-plot(post_traj, plot.type = "medianHDI")
-legend("topleft", legend=c("posterior median", "95% HDI"), lty=1, col=c("red","grey"))
+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')
 
 ```
-Alternatively we can plot an ensemble of posterior trajectories
+Again, this confirms that the inference procedure does indeed recover the true model.
+
+
+Lastly, we can plot an ensemble of posterior trajectories using `plot.type="ensemble"` 
 ```{r post-sims-ensemble}
 plot(post_traj, plot.type = "ensemble", col = "#FF000040")
 ```
+
+# References
diff --git a/inst/doc/logistic_ode_example.pdf b/inst/doc/logistic_ode_example.pdf
index 4cc33ee..03f3327 100644
Binary files a/inst/doc/logistic_ode_example.pdf and b/inst/doc/logistic_ode_example.pdf differ
diff --git a/vignettes/chytrid_dede_example.R b/vignettes/chytrid_dede_example.R
deleted file mode 100644
index c8c8986..0000000
--- a/vignettes/chytrid_dede_example.R
+++ /dev/null
@@ -1,165 +0,0 @@
-## ----install1, eval=FALSE------------------------------------------------
-#  install.packages("devtools")
-
-## ----install1b-----------------------------------------------------------
-#Load the devtools package.
-library(devtools)
-
-## ----install2, eval=FALSE, results='hide', message=FALSE, warning=FALSE----
-#  install_github("pboesu/debinfer")
-
-## ----loadlib, message=FALSE----------------------------------------------
-library(deBInfer)
-
-## ----dde-def-------------------------------------------------------------
-#dede version
-CSZ.dede<-function(t,y,p){
-
-  sr    <-p["sr"]
-  fs    <-p["fs"]
-  ds    <-p["ds"]
-  eta   <-p["eta"]
-  Tmin  <-p["Tmin"]
-  ##Tmax  <-p["Tmax"]
-  muz <-p["muz"]
-
-  Rs<-Ms<-0
-  lag1<-lag2<-0
-
-  if (t>Tmin){
-    lag1<-lagvalue(t-Tmin)
-    Rs <- sr*fs*lag1[1]
-  }
-
-  phiZ <- eta*y[2]
-  dy1 <- -(muz+sr)*y[1]
-  dy2 <- Rs - Ms - ds*y[2]
-  dy3 <- phiZ - (muz+sr)*y[3]
-
-  if(y[1]<0) dy1<-0
-  if(y[2]<0){
-    dy2 <- Rs - Ms
-    dy3 <- -(muz+sr)*y[3]
-  }
-  if(y[3]<0){
-    dy3 <- dy3+(muz+sr)*y[3]
-  }
-
-  list(c(dy1,dy2,dy3))
-}
-
-## ----data----------------------------------------------------------------
-#load chytrid data
-data(chytrid)
-#have a look at the variables
-head(chytrid)
-#plot the data
-plot(chytrid, xlab='Time (days)', ylab='Zoospores x 10e4', xlim=c(0,10))
-
-## ----obs-model-----------------------------------------------------------
-# observation model
-chytrid_obs_model<-function(data, sim.data, samp){
-
-  ec<-0.01
-  llik.Z<-0
-  for(i in unique(data$time)){
-    try(llik.Z<-llik.Z + sum(dpois(data$count[data$time==i], 
-                                   lambda=(sim.data[,'Z'][sim.data[,'time']==i]+ec),
-                                   log=TRUE)))
-  }
-  llik<-llik.Z
-  return(llik)
-}
-
-## ----vars----------------------------------------------------------------
-sr <- debinfer_par(name = "sr", var.type = "de", fixed = FALSE,
-                   value = 2, prior="gamma", hypers=list(shape = 5, rate = 1),
-                   prop.var=c(3,4), samp.type="rw-unif")
-
-fs <- debinfer_par(name = "fs", var.type = "de", fixed = FALSE,
-                   value = 0.5, prior="beta", hypers=list(shape1 = 1, shape2 = 1),
-                   prop.var=0.01, samp.type="ind")
-
-ds <- debinfer_par(name = "ds", var.type = "de", fixed = FALSE,
-                   value = 2, prior="gamma", hypers=list(shape = 1, rate = 1),
-                   prop.var=0.1, samp.type="rw")
-
-muz <- debinfer_par(name = "muz", var.type = "de", fixed = FALSE,
-                    value = 1, prior="gamma", hypers=list(shape = 5, rate = 1),
-                    prop.var=c(4,5), samp.type="rw-unif")
-
-eta <- debinfer_par(name = "eta", var.type = "de", fixed = FALSE,
-                    value = 10, prior="gamma", hypers=list(shape = 1, rate = 0.25),
-                    prop.var=5, samp.type="rw")
-
-Tmin <- debinfer_par(name = "Tmin", var.type = "de", fixed = FALSE,
-                     value = 3, prior="unif", hypers=list(min = 2, max = 6),
-                     prop.var=0.05, samp.type="rw")
-
-
-# ----inits---------------------------------------------------------------
-C <- debinfer_par(name = "C", var.type = "init", fixed = TRUE, value = 120)
-S <- debinfer_par(name = "S", var.type = "init", fixed = TRUE, value = 0)
-Z <- debinfer_par(name = "Z", var.type = "init", fixed = TRUE, value = 0)
-
-## ----mcmc-setup----------------------------------------------------------
-# ----setup---------------------------------------------------------------
-mcmc.pars <- setup_debinfer(sr, fs, ds, muz, eta, Tmin, C, S, Z)
-
-## ----de_mcmc, results='hide'---------------------------------------------
-# do inference with deBInfer
-# MCMC iterations
-iter = 500
-# inference call
-dede_rev <- de_mcmc(N = iter, data=chytrid, de.model=CSZ.dede,
-                               obs.model=chytrid_obs_model, all.params=mcmc.pars,
-                               Tmax = max(chytrid$time), data.times=c(0,chytrid$time), cnt=50,
-                               plot=FALSE, sizestep=0.1, solver="dede", verbose = TRUE)
-
-## ----message=FALSE, warning=FALSE,fig.width = 8, fig.height = 8----------
-par(mfrow = c(3,4))
-plot(dede_rev, ask=FALSE, auto.layout=FALSE)
-
-## ---- fig.width = 8, fig.height = 8--------------------------------------
-burnin = 50
-pairs(dede_rev, burnin = burnin, scatter=TRUE, trend=TRUE)
-post_prior_densplot(dede_rev, burnin = burnin)
-
-## ------------------------------------------------------------------------
-par(mfrow=c(2,3), mgp=c(2.2, 0.8, 0))
-#define a fancy y axis label
-ylabel = expression(paste(Pr,"(", theta,"|", "Y", ")"))
-#plot the individual parameters
-post_prior_densplot(dede_rev, param="sr",xlab=expression(theta), 
-                    ylab=ylabel, show.obs=FALSE, xlim=c(0,8), 
-                    main=expression(paste("s",phantom()[{paste("r")}])))
-legend("topright", legend=c("Posterior","Prior"), lty = 1, col = c("black", "red"))
-post_prior_densplot(dede_rev, param="fs",xlab=expression(theta), 
-                    ylab=ylabel, show.obs=FALSE, xlim=c(-0.1,1.1), 
-                    main=expression(paste("f",phantom()[{paste("s")}])))
-post_prior_densplot(dede_rev, param="ds",xlab=expression(theta),
-                    ylab=ylabel, show.obs=FALSE, xlim=c(0,3), 
-                    main=expression(paste("d",phantom()[{paste("s")}])))
-post_prior_densplot(dede_rev, param="muz",xlab=expression(theta),
-                    ylab=ylabel, show.obs=FALSE, xlim=c(0,6), 
-                    main=expression(paste(mu,phantom()[{paste("Z")}])))
-post_prior_densplot(dede_rev, param="eta",xlab=expression(theta),
-                    ylab=ylabel, show.obs=FALSE, xlim=c(0,50), ylim=c(0,0.2), 
-                    main=expression(eta))
-post_prior_densplot(dede_rev, param="Tmin",xlab=expression(theta),
-                    ylab=ylabel, show.obs=FALSE, xlim=c(1.5,6.5), 
-                    main=expression(paste("T",phantom()[{paste("min")}])))
-
-## ----post-sims-----------------------------------------------------------
-post_traj <- post_sim(dede_rev, n=100, times=0:100, burnin=burnin, output = 'all', prob = 0.95)
-
-## ----post-sims-plot, fig.width = 8, fig.height = 3-----------------------
-#median and HDI
-par(mfrow=c(1,3))
-plot(post_traj, plot.type = "medianHDI", auto.layout = FALSE)
-legend("topright", legend=c("posterior median", "95% HDI"), lty=1, col=c("red","grey"), bty='n')
-
-
-## ----post-sims-ensemble, fig.width = 8, fig.height = 6-------------------
-plot(post_traj, plot.type = "ensemble", col = "#FF000040")
-
diff --git a/vignettes/chytrid_dede_example.Rmd b/vignettes/chytrid_dede_example.Rmd
index 092e35d..1f78341 100644
--- a/vignettes/chytrid_dede_example.Rmd
+++ b/vignettes/chytrid_dede_example.Rmd
@@ -194,7 +194,7 @@ plot(dede_rev, ask=FALSE, auto.layout=FALSE)
 ```
 From the traceplot we can see that with only 500 iterations the chains have neither mixed well, nor reached stationarity. For demonstration purposes we remove a burnin period of 100 samples and have a look at parameter correlations, and the overlap between the posterior and prior densities.
 ```{r, fig.width = 8, fig.height = 8}
-burnin = 50
+burnin = 100
 pairs(dede_rev, burnin = burnin, scatter=TRUE, trend=TRUE)
 post_prior_densplot(dede_rev, burnin = burnin)
 ```
diff --git a/vignettes/chytrid_dede_example.pdf b/vignettes/chytrid_dede_example.pdf
deleted file mode 100644
index b142907..0000000
Binary files a/vignettes/chytrid_dede_example.pdf and /dev/null differ
diff --git a/vignettes/logistic_ode_example.R b/vignettes/logistic_ode_example.R
deleted file mode 100644
index a2b125e..0000000
--- a/vignettes/logistic_ode_example.R
+++ /dev/null
@@ -1,118 +0,0 @@
-## ----install1, eval=FALSE------------------------------------------------
-#  install.packages("devtools")
-
-## ----install1b-----------------------------------------------------------
-#Load the devtools package.
-library(devtools)
-
-## ----install2, eval=FALSE, results='hide', message=FALSE, warning=FALSE----
-#  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"--------------------------------------------
-# 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, 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)
-
-## ------------------------------------------------------------------------
-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")
-
diff --git a/vignettes/logistic_ode_example.Rmd b/vignettes/logistic_ode_example.Rmd
index 9da1286..0b6eb7b 100644
--- a/vignettes/logistic_ode_example.Rmd
+++ b/vignettes/logistic_ode_example.Rmd
@@ -172,7 +172,7 @@ 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)
 ```
-Overlay the true underlying parameter values as vertical dashed lines shows that our estimation is doing a pretty good job.
+Overlay the true underlying parameter values as vertical dashed lines shows that our estimation is recovering the true parameters.
 
 # Posterior simulations
 We then simulate DE model trajectories from the posterior and calculate the 95% highest posterior density interval for the deterministic part of the model.
diff --git a/vignettes/logistic_ode_example.pdf b/vignettes/logistic_ode_example.pdf
deleted file mode 100644
index 092d34d..0000000
Binary files a/vignettes/logistic_ode_example.pdf and /dev/null differ