## This example shows how to run and calibrate the VSEM model
library(BayesianTools)
# Create input data for the model
PAR <- VSEMcreatePAR(1:1000)
plotTimeSeries(observed = PAR)
# load reference parameter definition (upper, lower prior)
refPars <- VSEMgetDefaults()
# this adds one additional parameter for the likelihood standard deviation (see below)
refPars[12,] <- c(0.1, 0.001, 0.5)
rownames(refPars)[12] <- "error-sd"
head(refPars)
# create some simulated test data
# generally recommended to start with simulated data before moving to real data
referenceData <- VSEM(refPars$best[1:11], PAR) # model predictions with reference parameters
referenceData[,1] = 1000 * referenceData[,1]
# this adds the error - needs to conform to the error definition in the likelihood
obs <- referenceData + rnorm(length(referenceData), sd = refPars$best[12])
oldpar <- par(mfrow = c(2,2))
for (i in 1:4) plotTimeSeries(observed = obs[,i],
predicted = referenceData[,i], main = colnames(referenceData)[i])
# Best to program in a way that we can choose easily which parameters to calibrate
parSel = c(1:6, 12)
# here is the likelihood
likelihood <- function(x, sum = TRUE){
# createMixWithDefaults sets the parameters that are not calibrated on default values
x <- createMixWithDefaults(x, refPars$best, parSel)
predicted <- VSEM(x[1:11], PAR) # replace here VSEM with your model
predicted[,1] = 1000 * predicted[,1] # this is just rescaling
diff <- c(predicted[,1:4] - obs[,1:4]) # difference betweeno observed and predicted
# univariate normal likelihood. Note that there is a parameter involved here that is fit
llValues <- dnorm(diff, sd = x[12], log = TRUE)
if (sum == FALSE) return(llValues)
else return(sum(llValues))
}
# optional, you can also directly provide lower, upper in the createBayesianSetup, see help
prior <- createUniformPrior(lower = refPars$lower[parSel],
upper = refPars$upper[parSel], best = refPars$best[parSel])
bayesianSetup <- createBayesianSetup(likelihood, prior, names = rownames(refPars)[parSel])
# settings for the sampler, iterations should be increased for real applicatoin
settings <- list(iterations = 2000, nrChains = 2)
out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "DEzs", settings = settings)
plot(out)
summary(out)
gelmanDiagnostics(out) # should be below 1.05 for all parameters to demonstrate convergence
Run the code above in your browser using DataLab