Creates a standardized collection of prior, likelihood and posterior functions, including error checks etc.

Creates a standardized collection of prior, likelihood and posterior functions, including error checks etc.

createBayesianSetup(likelihood, prior = NULL, priorSampler = NULL,
  parallel = FALSE, lower = NULL, upper = NULL, best = NULL,
  names = NULL, parallelOptions = list(variables = "all", packages = "all",
  dlls = NULL), catchDuplicates = FALSE)

Arguments

likelihood	log likelihood density function
prior	either a prior class (see createPrior) or a log prior density
priorSampler	if a prior density (and not a prior class) is provided to prior, the optional prior sampling function can be provided here
parallel	parallelization option. Default is F. Other options are T, or "external". See details.
lower	vector with lower prior limits
upper	vector with upper prior limits
best	vector with best values
names	optional vector with parameter names
parallelOptions	list containing three lists. First "packages" determines the R packages necessary to run the likelihood function. Second "variables" the objects in the global envirnment needed to run the likelihood function and third "dlls" the DLLs needed to run the likelihood function (see Details and Examples).
catchDuplicates	Logical, determines whether unique parameter combinations should only be evaluated once. Only used when the likelihood accepts a matrix with parameter as columns.

Details

If lower and upper (and optionally best) but not a prior object are passed, the function will create an uniform prior with the limits lower and upper. If a prior object and lower and upper are passed, the function will ignore lower and upper, and will only use the prior object.

For parallelization, option T means that an automatic parallelization via R is attempted, or "external", in which case it is assumed that the likelihood is already parallelized. In this case it needs to accept a matrix with parameters as columns. Further you can specify the packages, objects and DLLs that are exported to the cluster. By default a copy of your workspace is exported. However, depending on your workspace this can be very inefficient.

For more details, make sure to read the vignette (run vignette("BayesianTools", package="BayesianTools")

See also

checkBayesianSetup createLikelihood createPrior

Examples

ll <- function(x) sum(dnorm(x, log = TRUE))

test <- createBayesianSetup(ll, prior = NULL, priorSampler = NULL, lower = -10, upper = 10)
str(test)
#> List of 7
#>  $ prior     :List of 6
#>   ..$ density        :function (x)  
#>   .. ..- attr(*, "srcref")=Class 'srcref'  atomic [1:8] 59 21 63 3 21 3 59 63
#>   .. .. .. ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' <environment: 0x000000002766d2b8> 
#>   ..$ sampler        :function (n = NULL)  
#>   .. ..- attr(*, "srcref")=Class 'srcref'  atomic [1:8] 68 24 76 5 24 5 68 76
#>   .. .. .. ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' <environment: 0x000000002766d2b8> 
#>   ..$ lower          : num -10
#>   ..$ upper          : num 10
#>   ..$ best           : num 0
#>   ..$ originalDensity:function (x)  
#>   .. ..- attr(*, "srcref")=Class 'srcref'  atomic [1:8] 101 14 104 3 14 3 101 104
#>   .. .. .. ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' <environment: 0x000000002766d2b8> 
#>   ..- attr(*, "class")= chr "prior"
#>  $ likelihood:List of 3
#>   ..$ density:function (x, ...)  
#>   .. ..- attr(*, "srcref")=Class 'srcref'  atomic [1:8] 45 21 95 3 21 3 45 95
#>   .. .. .. ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' <environment: 0x00000000255f6b60> 
#>   ..$ sampler: NULL
#>   ..$ cl     : NULL
#>   ..- attr(*, "class")= chr "likelihood"
#>  $ posterior :List of 1
#>   ..$ density:function (x, returnAll = F)  
#>   .. ..- attr(*, "srcref")=Class 'srcref'  atomic [1:8] 9 16 41 3 16 3 9 41
#>   .. .. .. ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' <environment: 0x0000000028eeb968> 
#>   ..- attr(*, "class")= chr "posterior"
#>  $ names     : chr "par 1"
#>  $ numPars   : int 1
#>  $ model     : NULL
#>  $ parallel  : logi FALSE
#>  - attr(*, "class")= chr "BayesianSetup"
test$prior$density(0)
#> [1] -2.995732

test$likelihood$density(c(1,1))
#> [1] -2.837877
test$likelihood$density(1)
#> [1] -1.418939
test$posterior$density(1)
#> [1] -4.414671
test$posterior$density(1, returnAll = TRUE)
#> [1] -4.414671 -1.418939 -2.995732

test$likelihood$density(matrix(rep(1,4), nrow = 2))
#> [1] -2.837877 -2.837877
#test$posterior$density(matrix(rep(1,4), nrow = 2), returnAll = TRUE)
test$likelihood$density(matrix(rep(1,4), nrow = 4))
#> [1] -1.418939 -1.418939 -1.418939 -1.418939



## Example of how to use parallelization using the VSEM model
# Note that the parallelization produces overhead and is not always
# speeding things up. In this example, due to the small
# computational cost of the VSEM the parallelization is
# most likely to reduce the speed of the sampler.

# Creating reference data
PAR <- VSEMcreatePAR(1:1000)
refPars   <- VSEMgetDefaults()
refPars[12,] <- c(0.2, 0.001, 1)
rownames(refPars)[12] <- "error-sd"

referenceData <- VSEM(refPars$best[1:11], PAR)
obs = apply(referenceData, 2, function(x) x + rnorm(length(x),
                                                    sd = abs(x) * refPars$best[12]))

# Selecting parameters
parSel = c(1:6, 12)


## Builidng the likelihood function
likelihood <- function(par, sum = TRUE){
  x = refPars$best
  x[parSel] = par
  predicted <- VSEM(x[1:11], PAR)
  diff = c(predicted[,1:3] - obs[,1:3])
  llValues = dnorm(diff, sd = max(abs(c(predicted[,1:3])),0.0001) * x[12], log = TRUE)
  if (sum == False) return(llValues)
  else return(sum(llValues))
}

# Prior
prior <- createUniformPrior(lower = refPars$lower[parSel], upper = refPars$upper[parSel])


####
## Definition of the packages and objects that are exported to the cluster.
# These are the objects that are used in the likelihood function.
opts <- list(packages = list("BayesianTools"), variables = list("refPars", "obs", "PAR" ),
             dlls = NULL)


# Create Bayesian Setup
BSVSEM <- createBayesianSetup(likelihood, prior, best = refPars$best[parSel],
                              names = rownames(refPars)[parSel], parallel = 2,
                              parallelOptions = opts)
#> Error in get(name, envir = envir): object 'refPars' not found

## The bayesianSetup can now be used in the runMCMC function.
# Note that not all samplers can make use of parallel
# computing.

# Remove the Bayesian Setup and close the cluster
stopParallel(BSVSEM)
rm(BSVSEM)
#> Warning: object 'BSVSEM' not found