# TODO change this once the mcmc branch is merged in! 
devtools::install_github('kdorheim/MENDplus', ref = 'mcmc_dev')
## Skipping install of 'MENDplus' from a github remote, the SHA1 (3d8384d4) has not changed since last install.
##   Use `force = TRUE` to force installation
library(MENDplus)
library(mcmc)

Objective:

In this vignette we walk through the code that is needed to set up the mcmc metrop function to to calibrate MENDplus.

A general over view of this process is the following

  1. Set up comparison data
  2. Deinfe the log posterior function
  3. Run the metropolis algorithm

1. Comparison Data

As an example we are going to use output from a MEND run as comparison data. Users can use data from observations or other modeling sources as comparison data. But it will need to be a data frame consisiting of a time and IC columns, where IC is the \(CO_[2]\).

# Specify the time steps.
times <- seq(0, 1000, 1)

# Define the inital state of the MEND C pools. 
state <- c(P = 10,  M = 5,  Q = 0.1,  B = 2,  D = 1,  EP = 0.00001, 
           EM = 0.00001, IC = 0,  Tot = 18.10002)


# Set up the ODE solver to solve MEND governed with MEND's default flux.
output <- as.data.frame(deSolve::ode(y = state, times = times, parms = default_parameters,   
                                     func = MEND_carbon_pools, flux_function = MEND_fluxes))  

comp <- as.data.frame(output[, names(output) %in% c('time', 'IC')])

2. Deinfe the log posterior function

The posterior probability to a key part of any Baeysian or MCMC calibration see here for more details. The metrop function requires

If a function, it evaluates the log unnormalized probability density of the desired equilibrium distribution of the Markov chain.

We devloped make_logpost that returns a function that calculates the log unnormalized probability based the user defined parameters to calibrate and comparison data.

# Select the MEND parameters to calibrate and provide an inital estimate for the values. 
params_to_calibrate <- c('K.p' = 25,'K.m' = 125, 'K.d' = 0.13)

# Use the make_logpost function to create the log posteior function.
log_post <- make_logpost(comp = comp, params = params_to_calibrate)

3. Run the mcmc!

## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 33.2342, R2 = 1.3829e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 33.2342, R2 = 1.3829e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 33.2342, R2 = 1.3829e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 33.2342, R2 = 1.3829e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 33.2342, R2 = 1.3829e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 33.2342, R2 = 1.3829e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 33.2342, R2 = 1.3829e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 33.2342, R2 = 1.3829e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 33.2342, R2 = 1.3829e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 33.2342, R2 = 1.3829e-15
##  
## DLSODA-  Above warning has been issued I1 times.  
##      It will not be issued again for this problem.
## In above message, I1 = 10
##  
## DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
##       taken on this call before reaching TOUT     
## In above message, I1 = 5000
##  
## In above message, R1 = 33.2342
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 20.4399, R2 = 5.36736e-16
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 20.4399, R2 = 5.36736e-16
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 20.4399, R2 = 5.36736e-16
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 20.4399, R2 = 5.36736e-16
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 20.4399, R2 = 1.07347e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 20.4399, R2 = 1.07347e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 20.4399, R2 = 1.07347e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 20.4399, R2 = 2.26499e-16
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 20.4399, R2 = 2.26499e-16
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 20.4399, R2 = 4.52997e-16
##  
## DLSODA-  Above warning has been issued I1 times.  
##      It will not be issued again for this problem.
## In above message, I1 = 10
##  
## DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
##       taken on this call before reaching TOUT     
## In above message, I1 = 5000
##  
## In above message, R1 = 20.4399
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 274.954, R2 = 1.9069e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 274.954, R2 = 1.9069e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 274.954, R2 = 1.9069e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 274.954, R2 = 1.9069e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 274.954, R2 = 2.00567e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 274.954, R2 = 2.00567e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 274.954, R2 = 2.00567e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 274.954, R2 = 2.00567e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 274.954, R2 = 2.00567e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 274.954, R2 = 2.00567e-14
##  
## DLSODA-  Above warning has been issued I1 times.  
##      It will not be issued again for this problem.
## In above message, I1 = 10
##  
## DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
##       taken on this call before reaching TOUT     
## In above message, I1 = 5000
##  
## In above message, R1 = 274.954
##  
## DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
##       taken on this call before reaching TOUT     
## In above message, I1 = 5000
##  
## In above message, R1 = 3.57631
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 421.594, R2 = 1.43086e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 421.594, R2 = 1.43086e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 421.594, R2 = 1.43086e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 421.594, R2 = 2.86172e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 421.594, R2 = 2.86172e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 421.594, R2 = 2.86172e-15
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 421.594, R2 = 1.83651e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 421.594, R2 = 1.83651e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 421.594, R2 = 1.83651e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 421.594, R2 = 1.32728e-14
##  
## DLSODA-  Above warning has been issued I1 times.  
##      It will not be issued again for this problem.
## In above message, I1 = 10
##  
## DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
##       taken on this call before reaching TOUT     
## In above message, I1 = 5000
##  
## In above message, R1 = 421.594
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 378.589, R2 = 2.41016e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 378.589, R2 = 2.41016e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 378.589, R2 = 2.41016e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 378.589, R2 = 2.41016e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 378.589, R2 = 2.41016e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 378.589, R2 = 2.41016e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 378.589, R2 = 2.41016e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 378.589, R2 = 2.41016e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 378.589, R2 = 2.41016e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 378.589, R2 = 2.41016e-14
##  
## DLSODA-  Above warning has been issued I1 times.  
##      It will not be issued again for this problem.
## In above message, I1 = 10
##  
## DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
##       taken on this call before reaching TOUT     
## In above message, I1 = 5000
##  
## In above message, R1 = 378.589
##  
## DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
##       taken on this call before reaching TOUT     
## In above message, I1 = 5000
##  
## In above message, R1 = 2.99161
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 547.026, R2 = 5.24296e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 547.026, R2 = 5.24296e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 547.026, R2 = 5.24296e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 547.026, R2 = 5.4415e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 547.026, R2 = 5.4415e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 547.026, R2 = 5.4415e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 547.026, R2 = 4.5081e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 547.026, R2 = 4.5081e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 547.026, R2 = 4.5081e-14
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 547.026, R2 = 1.08346e-14
##  
## DLSODA-  Above warning has been issued I1 times.  
##      It will not be issued again for this problem.
## In above message, I1 = 10
##  
## DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
##       taken on this call before reaching TOUT     
## In above message, I1 = 5000
##  
## In above message, R1 = 547.026
##  
## DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
##       taken on this call before reaching TOUT     
## In above message, I1 = 5000
##  
## In above message, R1 = 2.04546
##

Take a look at some quick diagnostic plots.

plot(ts(out$batch))

acf(out$batch)