# inst/examples/plotTimeSeriesResultsHelp.R In BayesianTools: General-Purpose MCMC and SMC Samplers and Tools for Bayesian Statistics

```# Create input data for the model
# see help for the VSEM model
PAR <- VSEMcreatePAR(1:1000)

# load reference parameter definition (upper, lower prior)
refPars <- VSEMgetDefaults()
# this adds one additional parameter for the likelihood standard deviation (see below)
refPars[12,] <- c(2, 0.1, 4)
rownames(refPars) <- "error-sd"

# 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)

parSel = c(1:6, 12)

# here is the likelihood
likelihood <- function(par, sum = TRUE){
# set parameters that are not calibrated on default values
x = refPars\$best
x[parSel] = par
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, 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)

\dontrun{

out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "DEzs", settings = settings)

# Posterior predictive simulations

# Create a prediction function
createPredictions <- function(par){
# set the parameters that are not calibrated on default values
x = refPars\$best
x[parSel] = par
predicted <- VSEM(x[1:11], PAR) # replace here VSEM with your model
return(predicted[,1] * 1000)
}

# Create an error function
createError <- function(mean, par){
return(rnorm(length(mean), mean = mean, sd = par))
}

# plot prior predictive distribution and prior predictive simulations
plotTimeSeriesResults(sampler = out, model = createPredictions, observed = obs[,1],
error = createError, prior = TRUE, main = "Prior predictive")

}
```

## Try the BayesianTools package in your browser

Any scripts or data that you put into this service are public.

BayesianTools documentation built on Dec. 10, 2019, 1:08 a.m.