## ----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 = 1000 settings = list(iterations = iter, message = FALSE) out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings) ## ----------------------------------------------------------------------------- out2 <- runMCMC(bayesianSetup = out) ## ----------------------------------------------------------------------------- opar = par(mfrow=c(2,1)) plot(getSample(out, whichParameters = 1), type = "l", ylab = "par1", xlab = "iteration") plot(getSample(out2, whichParameters = 1) , type = "l", ylab = "par1", xlab = "iteration") par(opar) ## ----echo = T----------------------------------------------------------------- settings = list(iterations = iter, nrChains = 3, message = FALSE) out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings) ## ----------------------------------------------------------------------------- tracePlot(out) # identical to plot(out) ## ----------------------------------------------------------------------------- gelmanDiagnostics(out, plot = T, start = 100) ## ----------------------------------------------------------------------------- print(out, start = 100) summary(out, start = 100) ## ----------------------------------------------------------------------------- correlationPlot(out) marginalPlot(out, prior = TRUE) ## ----------------------------------------------------------------------------- marginalLikelihood(out) DIC(out) MAP(out) ## ----eval = F----------------------------------------------------------------- # getSample(out, start = 100, end = NULL, thin = 5, whichParameters = 1:2) ## ----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 values # } # # ## Create Bayesian Setup # BS <- createBayesianSetup(likelihood, parallel = "external", ...) # # ## Run MCMC # runMCMC(BS, sampler = "SMC", ...) ## ----eval = FALSE------------------------------------------------------------- # ## n = Number of cores # n = 2 # x <- c(1:10) # likelihood <- function(param) return(sum(dnorm(x, mean = param, log = T))) # bayesianSetup <- createBayesianSetup(likelihood, parallel = n, lower = -5, upper = 5) # # ## give runMCMC a matrix with n rows of proposals as startValues or sample n times from the previous created sampler # out <- runMCMC(bayesianSetup, settings = list(iterations = 1000)) ## ----eval = FALSE------------------------------------------------------------- # ### Create cluster with n cores # cl <- parallel::makeCluster(n) # # ## Definition of the likelihood # likelihood <- function(X) sum(dnorm(c(1:10), mean = X, log = T)) # # ## Definition of the likelihood which will be calculated in parallel. Instead of the parApply function, we could also define a costly parallelized likelihood # pLikelihood <- function(param) parallel::parApply(cl = cl, X = param, MARGIN = 1, FUN = likelihood) # # ## export functions, dlls, libraries # # parallel::clusterEvalQ(cl, library(BayesianTools)) # parallel::clusterExport(cl, varlist = list(likelihood)) # # ## create BayesianSetup # bayesianSetup <- createBayesianSetup(pLikelihood, lower = -10, upper = 10, parallel = 'external') # # ## For this case we want to parallelize the internal chains, therefore we create a n row matrix with startValues, if you parallelize a model in the likelihood, do not set a n*row Matrix for startValue # settings = list(iterations = 100, nrChains = 1, startValue = bayesianSetup$prior$sampler(n)) # # ## runMCMC # out <- runMCMC(bayesianSetup, settings, sampler = "DEzs") ## ----eval = FALSE------------------------------------------------------------- # ### Create cluster with n cores # cl <- parallel::makeCluster(n) # # ## export your model # # parallel::clusterExport(cl, varlist = list(complexModel)) # # ## Definition of the likelihood # likelihood <- function(param) { # # ll <- complexModel(param) # # return(ll) # } # # ## create BayesianSetup and settings # bayesianSetup <- createBayesianSetup(likelihood, lower = -10, upper = 10, parallel = 'external') # settings = list(iterations = 100, nrChains = 1) # # ## runMCMC # out <- runMCMC(bayesianSetup, settings) # ## ----eval = FALSE------------------------------------------------------------- # ### Definition of likelihood function # x <- c(1:10) # likelihood <- function(param) return(sum(dnorm(x, mean = param, log = T))) # # ## Create BayesianSetup and settings # bayesianSetup <- createBayesianSetup(likelihood, lower = -10, upper = 10, parallel = F) # settings = list(iterations = 100000) # # ## Start cluster with n cores for n chains and export BayesianTools library # cl <- parallel::makeCluster(n) # parallel::clusterEvalQ(cl, library(BayesianTools)) # # ## calculate parallel n chains, for each chain the likelihood will be calculated on one core # out <- parallel::parLapply(cl, 1:n, fun = function(X, bayesianSetup, settings) runMCMC(bayesianSetup, settings, sampler = "DEzs"), bayesianSetup, settings) # # ## Combine the chains # out <- createMcmcSamplerList(out) ## ----eval = F----------------------------------------------------------------- # createBetaPrior(a, b, lower, upper) # createUniformPrior(lower, upper, best) # createTruncatedNormalPrior(mean, sd, lower, upper) ## ----------------------------------------------------------------------------- # 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) ## ----------------------------------------------------------------------------- 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) ## ----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.ML - M2$ln.ML) ## ----------------------------------------------------------------------------- exp(M1$ln.ML) / ( exp(M1$ln.ML) + exp(M2$ln.ML)) ## ----------------------------------------------------------------------------- DIC(out1)$DIC DIC(out2)$DIC ## ----------------------------------------------------------------------------- # This will not work, since likelihood1 has no sum argument # WAIC(out1) # likelihood with sum argument likelihood3 <- function(param, sum = TRUE){ pred <- param[1] + param[2]*x + param[3] * x^2 singlelikelihoods <- dnorm(y, mean = pred, sd = 1/(param[4]^2), log = T) return(if (sum == TRUE) sum(singlelikelihoods) else singlelikelihoods) } setUp3 <- createBayesianSetup(likelihood3, lower = c(-5,-5,-5,0.01), upper = c(5,5,5,30)) out3 <- runMCMC(bayesianSetup = setUp3, sampler = "Metropolis", settings = settings) WAIC(out3)