## ----global_options, include=FALSE---------------------------------------
knitr::opts_chunk$set(fig.width=5, fig.height=5, warning=FALSE, cache = F)

## ---- echo = F, message = F----------------------------------------------
set.seed(123)

## ---- eval = F-----------------------------------------------------------
#  install.packages("BayesianTools")

## ------------------------------------------------------------------------
library(BayesianTools)
citation("BayesianTools")

## ------------------------------------------------------------------------
set.seed(123)

## ---- eval = F-----------------------------------------------------------
#  sessionInfo()

## ------------------------------------------------------------------------
ll <- generateTestDensityMultiNormal(sigma = "no correlation")
bayesianSetup = createBayesianSetup(likelihood = ll, lower = rep(-10, 3), upper = rep(10, 3))

## ------------------------------------------------------------------------
iter = 10000
settings = list(iterations = iter, message = FALSE)
out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)

## ------------------------------------------------------------------------
print(out)
summary(out)

## ------------------------------------------------------------------------
plot(out) # plot internally calls tracePlot(out)
correlationPlot(out)
marginalPlot(out)

## ------------------------------------------------------------------------
marginalLikelihood(out)
DIC(out)
MAP(out)

## ---- eval = F-----------------------------------------------------------
#  getSample(out, start = 100, end = NULL, thin = 5, whichParameters = 1:2)

## ---- echo = T-----------------------------------------------------------
iter = 1000
settings = list(iterations = iter, nrChains = 3, message = FALSE)
out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)


## ------------------------------------------------------------------------
print(out)
summary(out)

## ------------------------------------------------------------------------
plot(out)

## ------------------------------------------------------------------------
#getSample(out, coda = F)
gelmanDiagnostics(out, plot = T)

## ---- eval = F-----------------------------------------------------------
#  ll = logDensity(x)

## ------------------------------------------------------------------------
ll <- generateTestDensityMultiNormal(sigma = "no correlation")
bayesianSetup = createBayesianSetup(likelihood = ll, lower = rep(-10, 3), upper = rep(10, 3))

## ---- eval = FALSE-------------------------------------------------------
#  ## Definition of likelihood function
#  likelihood <- function(matrix){
#  	# Calculate likelihood in parallel
#  	# Return vector of likelihood valus
#  }
#  
#  ## Create Bayesian Setup
#  BS <- createBayesianSetup(likelihood, parallel = "external" ...)
#  
#  ## Run MCMC
#  runMCMC(BS, sampler = "SMC", ...)

## ------------------------------------------------------------------------
# Create a BayesianSetup
ll <- generateTestDensityMultiNormal(sigma = "no correlation")
bayesianSetup = createBayesianSetup(likelihood = ll, lower = rep(-10, 3), upper = rep(10, 3))

settings = list(iterations = 2500,  message = FALSE)
out <- runMCMC(bayesianSetup = bayesianSetup, settings = settings)


newPrior = createPriorDensity(out, method = "multivariate", eps = 1e-10, lower = rep(-10, 3), upper =  rep(10, 3), best = NULL)
bayesianSetup <- createBayesianSetup(likelihood = ll, prior = newPrior)

settings = list(iterations = 1000,  message = FALSE)
out <- runMCMC(bayesianSetup = bayesianSetup, settings = settings)

## ---- message = F--------------------------------------------------------

ll <- generateTestDensityMultiNormal(sigma = "no correlation")

bayesianSetup = createBayesianSetup(likelihood = ll, lower = rep(-10, 3), upper = rep(10, 3))

settings = list(iterations = 10000, nrChains= 3, message = FALSE)
out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)

plot(out)
marginalPlot(out)
correlationPlot(out)
gelmanDiagnostics(out, plot=T)

# option to restart the sampler

settings = list(iterations = 1000, nrChains= 1, message = FALSE)
out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)

out2 <- runMCMC(bayesianSetup = out)

out3 <- runMCMC(bayesianSetup = out2)

#plot(out)
#plot(out3)

# create new prior from posterior sample 

newPriorFromPosterior <- createPriorDensity(out2)


## ------------------------------------------------------------------------
iter = 10000

## ------------------------------------------------------------------------
applySettingsDefault(sampler = "Metropolis")

## ---- results = 'hide', eval = F-----------------------------------------
#  settings <- list(iterations = iter, adapt = F, DRlevels = 1, gibbsProbabilities = NULL, temperingFunction = NULL, optimize = F, message = FALSE)
#  out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
#  plot(out)

## ---- results = 'hide', eval = F-----------------------------------------
#  settings <- list(iterations = iter, adapt = F, DRlevels = 1, gibbsProbabilities = NULL, temperingFunction = NULL, optimize = T, message = FALSE)
#  out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
#  plot(out)

## ---- results = 'hide', eval = F-----------------------------------------
#  settings <- list(iterations = iter, adapt = T, DRlevels = 1, gibbsProbabilities = NULL, temperingFunction = NULL, optimize = T,  message = FALSE)
#  out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
#  plot(out)

## ---- results = 'hide', eval = F-----------------------------------------
#  settings <- list(iterations = iter, adapt = F, DRlevels = 2, gibbsProbabilities = NULL, temperingFunction = NULL, optimize = T,  message = FALSE)
#  out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
#  plot(out)

## ---- results = 'hide', eval = F-----------------------------------------
#  settings <- list(iterations = iter, adapt = T, DRlevels = 2, gibbsProbabilities = NULL, temperingFunction = NULL, optimize = T,  message = FALSE)
#  out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
#  plot(out)

## ---- results = 'hide', eval = F-----------------------------------------
#  settings <- list(iterations = iter, adapt = T, DRlevels = 1, gibbsProbabilities = c(1,0.5,0), temperingFunction = NULL, optimize = T,  message = FALSE)
#  out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
#  plot(out)

## ---- results = 'hide', eval = F-----------------------------------------
#  temperingFunction <- function(x) 5 * exp(-0.01*x) + 1
#  settings <- list(iterations = iter, adapt = F, DRlevels = 1, gibbsProbabilities = c(1,1,0), temperingFunction = temperingFunction, optimize = T,  message = FALSE)
#  out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
#  plot(out)

## ---- results = 'hide', eval = F-----------------------------------------
#  settings <- list(iterations = iter,  message = FALSE)
#  out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "DE", settings = settings)
#  plot(out)

## ---- results = 'hide', eval = F-----------------------------------------
#  settings <- list(iterations = iter,  message = FALSE)
#  out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "DEzs", settings = settings)
#  plot(out)

## ---- results = 'hide', eval = F-----------------------------------------
#  settings <- list(iterations = iter,  message = FALSE)
#  out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "DREAM", settings = settings)
#  plot(out)

## ---- results = 'hide', eval = FALSE-------------------------------------
#  settings <- list(iterations = iter,  message = FALSE)
#  out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "DREAMzs", settings = settings)
#  #plot(out)

## ---- eval = F-----------------------------------------------------------
#  settings = list(iterations = iter,  message = FALSE)
#  
#  out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Twalk", settings = settings)

## ---- eval = T-----------------------------------------------------------
settings <- list(iterations = iter, nrChains = 3,  message = FALSE)
out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings)
plot(out)

#chain = getSample(out, coda = T)
gelmanDiagnostics(out, plot = F)

## ---- results = 'hide', eval = F-----------------------------------------
#  settings <- list(initialParticles = iter, iterations= 1)
#  out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "SMC", settings = settings)
#  plot(out)

## ---- results = 'hide', eval = F-----------------------------------------
#  settings <- list(initialParticles = iter, iterations= 10)
#  out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "SMC", settings = settings)
#  plot(out)

## ------------------------------------------------------------------------
sampleSize = 30
x <- (-(sampleSize-1)/2):((sampleSize-1)/2)
y <-  1 * x + 1*x^2 + rnorm(n=sampleSize,mean=0,sd=10)
plot(x,y, main="Test Data")

## ------------------------------------------------------------------------
likelihood1 <- function(param){
    pred = param[1] + param[2]*x + param[3] * x^2
    singlelikelihoods = dnorm(y, mean = pred, sd = 1/(param[4]^2), log = T)
    return(sum(singlelikelihoods))  
}

likelihood2 <- function(param){
    pred = param[1] + param[2]*x 
    singlelikelihoods = dnorm(y, mean = pred, sd = 1/(param[3]^2), log = T)
    return(sum(singlelikelihoods))  
}

## ------------------------------------------------------------------------
setUp1 <- createBayesianSetup(likelihood1, lower = c(-5,-5,-5,0.01), upper = c(5,5,5,30))

setUp2 <- createBayesianSetup(likelihood2, lower = c(-5,-5,0.01), upper = c(5,5,30))

## ---- results = 'hide'---------------------------------------------------
settings = list(iterations = 15000,  message = FALSE)
out1 <- runMCMC(bayesianSetup = setUp1, sampler = "Metropolis", settings = settings)
#tracePlot(out1, start = 5000)
M1 = marginalLikelihood(out1)
M1

settings = list(iterations = 15000,  message = FALSE)
out2 <- runMCMC(bayesianSetup = setUp2, sampler = "Metropolis", settings = settings)
#tracePlot(out2, start = 5000)
M2 = marginalLikelihood(out2)
M2

## ------------------------------------------------------------------------
exp(M1$ln.m - M2$ln.m)

## ------------------------------------------------------------------------
DIC(out1)$DIC
DIC(out2)$DIC

