inst/doc/bridgesampling_example_nimble.R

In bridgesampling: Bridge Sampling for Marginal Likelihoods and Bayes Factors

## -----------------------------------------------------------------------------
library(bridgesampling)

### generate data ###
set.seed(12345)

mu <- 0
tau2 <- 0.5
sigma2 <- 1

n <- 20
theta <- rnorm(n, mu, sqrt(tau2))
y <- rnorm(n, theta, sqrt(sigma2))
  

## ----eval=FALSE---------------------------------------------------------------
#  ### set prior parameters ###
#  mu0 <- 0
#  tau20 <- 1
#  alpha <- 1
#  beta <- 1

## ---- eval=FALSE--------------------------------------------------------------
#  library("nimble")
#  
#  # models
#  codeH0 <- nimbleCode({
#    invTau2 ~ dgamma(1, 1)
#    tau2 <- 1/invTau2
#    for (i in 1:20) {
#      theta[i] ~ dnorm(0, sd = sqrt(tau2))
#      y[i] ~ dnorm(theta[i], sd = 1)
#    }
#  })
#  codeH1 <- nimbleCode({
#    mu ~ dnorm(0, sd = 1)
#    invTau2 ~ dgamma(1, 1)
#    tau2 <- 1/invTau2
#    for (i in 1:20) {
#      theta[i] ~ dnorm(mu, sd = sqrt(tau2))
#      y[i] ~ dnorm(theta[i], sd = 1)
#    }
#  })
#  
#  ## steps for H0:
#  modelH0 <- nimbleModel(codeH0)
#  modelH0$setData(y = y) # set data
#  cmodelH0 <- compileNimble(modelH0) # make compiled version from generated C++
#  
#  ## steps for H1:
#  modelH1 <- nimbleModel(codeH1)
#  modelH1$setData(y = y) # set data
#  cmodelH1 <- compileNimble(modelH1) # make compiled version from generated C++
#  

## ---- eval=FALSE--------------------------------------------------------------
#  
#  # build MCMC functions, skipping customization of the configuration.
#  mcmcH0 <- buildMCMC(modelH0,
#                      monitors = modelH0$getNodeNames(stochOnly = TRUE,
#                                                      includeData = FALSE))
#  mcmcH1 <- buildMCMC(modelH1,
#                      monitors = modelH1$getNodeNames(stochOnly = TRUE,
#                                                      includeData = FALSE))
#  # compile the MCMC function via generated C++
#  cmcmcH0 <- compileNimble(mcmcH0, project = modelH0)
#  cmcmcH1 <- compileNimble(mcmcH1, project = modelH1)
#  
#  # run the MCMC.  This is a wrapper for cmcmc$run() and extraction of samples.
#  # the object samplesH1 is actually not needed as the samples are also in cmcmcH1
#  samplesH0 <- runMCMC(cmcmcH0, niter = 1e5, nburnin = 1000, nchains = 2,
#                       progressBar = FALSE)
#  samplesH1 <- runMCMC(cmcmcH1, niter = 1e5, nburnin = 1000, nchains = 2,
#                       progressBar = FALSE)

## ---- echo=FALSE--------------------------------------------------------------
load(system.file("extdata/", "vignette_example_nimble.RData",
                     package = "bridgesampling"))

## ----eval=FALSE---------------------------------------------------------------
#  # compute log marginal likelihood via bridge sampling for H0
#  H0.bridge <- bridge_sampler(cmcmcH0, silent = TRUE)
#  
#  # compute log marginal likelihood via bridge sampling for H1
#  H1.bridge <- bridge_sampler(cmcmcH1, silent = TRUE)

## -----------------------------------------------------------------------------
print(H0.bridge)
print(H1.bridge)

## ----eval=FALSE---------------------------------------------------------------
#  # compute percentage errors
#  H0.error <- error_measures(H0.bridge)$percentage
#  H1.error <- error_measures(H1.bridge)$percentage

## -----------------------------------------------------------------------------
print(H0.error)
print(H1.error)

## -----------------------------------------------------------------------------
# compute Bayes factor
BF01 <- bf(H0.bridge, H1.bridge)
print(BF01)

## -----------------------------------------------------------------------------
# compute posterior model probabilities (assuming equal prior model probabilities)
post1 <- post_prob(H0.bridge, H1.bridge)
print(post1)

## -----------------------------------------------------------------------------
# compute posterior model probabilities (using user-specified prior model probabilities)
post2 <- post_prob(H0.bridge, H1.bridge, prior_prob = c(.6, .4))
print(post2)