For today’s lab, we’ll use the BayesianTools
package, which has a bunch of different methods built in.
set.seed(7)
library(BayesianTools)
library(deSolve)
library(ggplot2)
Methods marked in italics are included in the BayesianTools
package, but I’ve included a couple of other common ones since they’re in commonly used packages like Stan
or POMP
. Links to more info included as well (also more info here)!
As a simple example to get to know the BayesianTools
package, let’s set up a similar normal distribution model to the one we built last time. We’ll suppose the data comes from a normal distribution with unknown \(\mu\), and known \(\sigma = 30\).
Let’s simulate some basic data to run our example:
data = rnorm(50,0,30) # here we're supposing the true mean is zero, and sampling 55 data points
The likelihood is just a normal distribution evaluated for our data set and whatever value of \(\mu\) we’re considering using:
loglike = function(par){
sum(dnorm(data, par, 30, log = T))
}
Now, we create a BayesianSetup—this is an object that has the prior and likelihood information that BayesianTools
needs.
We can set a prior using a distribution, or if we want to assume a uniform distribution with a given range, we can just provide the upper and lower limits. Note that if we don’t specify a prior, BayesianTools
will use an infinite uniform prior (basically equivalent to just using the likelihood).
We’ll also set the total number of iterations to run. We can set the starting values for the parameters, but if we don’t, BayesianTools
will default to sampling them from the prior, which is often a good idea anyhow.
setup = createBayesianSetup(likelihood = loglike, lower = -20, upper = 20)
settings = list(iterations = 10000, message=F)
Now, let’s run the MCMC! We’ll use the default algorithm for now (a form of differential evolution). Let’s also spit out some summary statistics and plots.
res = runMCMC(bayesianSetup = setup, settings = settings)
summary(res)
## # # # # # # # # # # # # # # # # # # # # # # # # #
## ## MCMC chain summary ##
## # # # # # # # # # # # # # # # # # # # # # # # # #
##
## # MCMC sampler: DEzs
## # Nr. Chains: 3
## # Iterations per chain: 3334
## # Rejection rate: 0.545
## # Effective sample size: 2295
## # Runtime: 0.891 sec.
##
## # Parameters
## psf MAP 2.5% median 97.5%
## par 1 1.006 7.158 -1.092 7.09 15.186
##
## ## DIC: 483.844
## ## Convergence
## Gelman Rubin multivariate psrf:
##
## ## Correlations
## par 1
## par 1 1
marginalPlot(res)
plot(res)