Skip to content

Commit

Permalink
rebuild vignettes
Browse files Browse the repository at this point in the history
  • Loading branch information
@singmann
singmann committed Aug 19, 2018
1 parent fb45daa commit af8990f
Show file tree
Hide file tree
Showing 9 changed files with 479 additions and 372 deletions.
264 changes: 137 additions & 127 deletions inst/doc/bridgesampling_example_jags.R
View file
Expand Up @@ -13,139 +13,149 @@ theta <- rnorm(n, mu, sqrt(tau2))
y <- rnorm(n, theta, sqrt(sigma2))


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

## ---- message = FALSE, results='hide'------------------------------------
library(R2jags)

### functions to get posterior samples ###

# H0: mu = 0
getSamplesModelH0 <- function(data, niter = 52000, nburnin = 2000, nchains = 3) {

model <- "
model {
for (i in 1:n) {
theta[i] ~ dnorm(0, invTau2)
y[i] ~ dnorm(theta[i], 1/sigma2)
}
invTau2 ~ dgamma(alpha, beta)
tau2 <- 1/invTau2
}"

s <- jags(data, parameters.to.save = c("theta", "invTau2"),
model.file = textConnection(model),
n.chains = nchains, n.iter = niter,
n.burnin = nburnin, n.thin = 1)

return(s)

}

# H1: mu != 0
getSamplesModelH1 <- function(data, niter = 52000, nburnin = 2000,
nchains = 3) {

model <- "
model {
for (i in 1:n) {
theta[i] ~ dnorm(mu, invTau2)
y[i] ~ dnorm(theta[i], 1/sigma2)
}
mu ~ dnorm(mu0, 1/tau20)
invTau2 ~ dgamma(alpha, beta)
tau2 <- 1/invTau2
}"

s <- jags(data, parameters.to.save = c("theta", "mu", "invTau2"),
model.file = textConnection(model),
n.chains = nchains, n.iter = niter,
n.burnin = nburnin, n.thin = 1)

return(s)

}

### get posterior samples ###

# create data lists for JAGS
data_H0 <- list(y = y, n = length(y), alpha = alpha, beta = beta, sigma2 = sigma2)
data_H1 <- list(y = y, n = length(y), mu0 = mu0, tau20 = tau20, alpha = alpha,
beta = beta, sigma2 = sigma2)

# fit models
samples_H0 <- getSamplesModelH0(data_H0)
samples_H1 <- getSamplesModelH1(data_H1)


## ------------------------------------------------------------------------
### functions for evaluating the unnormalized posteriors on log scale ###

log_posterior_H0 <- function(samples.row, data) {

mu <- 0
invTau2 <- samples.row[[ "invTau2" ]]
theta <- samples.row[ paste0("theta[", seq_along(data$y), "]") ]

sum(dnorm(data$y, theta, data$sigma2, log = TRUE)) +
sum(dnorm(theta, mu, 1/sqrt(invTau2), log = TRUE)) +
dgamma(invTau2, data$alpha, data$beta, log = TRUE)

}

log_posterior_H1 <- function(samples.row, data) {

mu <- samples.row[[ "mu" ]]
invTau2 <- samples.row[[ "invTau2" ]]
theta <- samples.row[ paste0("theta[", seq_along(data$y), "]") ]

sum(dnorm(data$y, theta, data$sigma2, log = TRUE)) +
sum(dnorm(theta, mu, 1/sqrt(invTau2), log = TRUE)) +
dnorm(mu, data$mu0, sqrt(data$tau20), log = TRUE) +
dgamma(invTau2, data$alpha, data$beta, log = TRUE)

}

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

## ---- eval=FALSE---------------------------------------------------------
# library(R2jags)
#
# ### functions to get posterior samples ###
#
# # H0: mu = 0
# getSamplesModelH0 <- function(data, niter = 52000, nburnin = 2000, nchains = 3) {
#
# model <- "
# model {
# for (i in 1:n) {
# theta[i] ~ dnorm(0, invTau2)
# y[i] ~ dnorm(theta[i], 1/sigma2)
# }
# invTau2 ~ dgamma(alpha, beta)
# tau2 <- 1/invTau2
# }"
#
# s <- jags(data, parameters.to.save = c("theta", "invTau2"),
# model.file = textConnection(model),
# n.chains = nchains, n.iter = niter,
# n.burnin = nburnin, n.thin = 1)
#
# return(s)
#
# }
#
# # H1: mu != 0
# getSamplesModelH1 <- function(data, niter = 52000, nburnin = 2000,
# nchains = 3) {
#
# model <- "
# model {
# for (i in 1:n) {
# theta[i] ~ dnorm(mu, invTau2)
# y[i] ~ dnorm(theta[i], 1/sigma2)
# }
# mu ~ dnorm(mu0, 1/tau20)
# invTau2 ~ dgamma(alpha, beta)
# tau2 <- 1/invTau2
# }"
#
# s <- jags(data, parameters.to.save = c("theta", "mu", "invTau2"),
# model.file = textConnection(model),
# n.chains = nchains, n.iter = niter,
# n.burnin = nburnin, n.thin = 1)
#
# return(s)
#
# }
#
# ### get posterior samples ###
#
# # create data lists for JAGS
# data_H0 <- list(y = y, n = length(y), alpha = alpha, beta = beta, sigma2 = sigma2)
# data_H1 <- list(y = y, n = length(y), mu0 = mu0, tau20 = tau20, alpha = alpha,
# beta = beta, sigma2 = sigma2)
#
# # fit models
# samples_H0 <- getSamplesModelH0(data_H0)
# samples_H1 <- getSamplesModelH1(data_H1)
#

## ----eval=FALSE----------------------------------------------------------
# ### functions for evaluating the unnormalized posteriors on log scale ###
#
# log_posterior_H0 <- function(samples.row, data) {
#
# mu <- 0
# invTau2 <- samples.row[[ "invTau2" ]]
# theta <- samples.row[ paste0("theta[", seq_along(data$y), "]") ]
#
# sum(dnorm(data$y, theta, data$sigma2, log = TRUE)) +
# sum(dnorm(theta, mu, 1/sqrt(invTau2), log = TRUE)) +
# dgamma(invTau2, data$alpha, data$beta, log = TRUE)
#
# }
#
# log_posterior_H1 <- function(samples.row, data) {
#
# mu <- samples.row[[ "mu" ]]
# invTau2 <- samples.row[[ "invTau2" ]]
# theta <- samples.row[ paste0("theta[", seq_along(data$y), "]") ]
#
# sum(dnorm(data$y, theta, data$sigma2, log = TRUE)) +
# sum(dnorm(theta, mu, 1/sqrt(invTau2), log = TRUE)) +
# dnorm(mu, data$mu0, sqrt(data$tau20), log = TRUE) +
# dgamma(invTau2, data$alpha, data$beta, log = TRUE)
#
# }
#

## ----eval=FALSE----------------------------------------------------------
# # specify parameter bounds H0
# cn <- colnames(samples_H0$BUGSoutput$sims.matrix)
# cn <- cn[cn != "deviance"]
# lb_H0 <- rep(-Inf, length(cn))
# ub_H0 <- rep(Inf, length(cn))
# names(lb_H0) <- names(ub_H0) <- cn
# lb_H0[[ "invTau2" ]] <- 0
#
# # specify parameter bounds H1
# cn <- colnames(samples_H1$BUGSoutput$sims.matrix)
# cn <- cn[cn != "deviance"]
# lb_H1 <- rep(-Inf, length(cn))
# ub_H1 <- rep(Inf, length(cn))
# names(lb_H1) <- names(ub_H1) <- cn
# lb_H1[[ "invTau2" ]] <- 0

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

## ----eval=FALSE----------------------------------------------------------
# # compute log marginal likelihood via bridge sampling for H0
# H0.bridge <- bridge_sampler(samples = samples_H0, data = data_H0,
# log_posterior = log_posterior_H0, lb = lb_H0,
# ub = ub_H0, silent = TRUE)
#
# # compute log marginal likelihood via bridge sampling for H1
# H1.bridge <- bridge_sampler(samples = samples_H1, data = data_H1,
# log_posterior = log_posterior_H1, lb = lb_H1,
# ub = ub_H1, silent = TRUE)

## ------------------------------------------------------------------------
# specify parameter bounds H0
cn <- colnames(samples_H0$BUGSoutput$sims.matrix)
cn <- cn[cn != "deviance"]
lb_H0 <- rep(-Inf, length(cn))
ub_H0 <- rep(Inf, length(cn))
names(lb_H0) <- names(ub_H0) <- cn
lb_H0[[ "invTau2" ]] <- 0

# specify parameter bounds H1
cn <- colnames(samples_H1$BUGSoutput$sims.matrix)
cn <- cn[cn != "deviance"]
lb_H1 <- rep(-Inf, length(cn))
ub_H1 <- rep(Inf, length(cn))
names(lb_H1) <- names(ub_H1) <- cn
lb_H1[[ "invTau2" ]] <- 0

## ------------------------------------------------------------------------
# compute log marginal likelihood via bridge sampling for H0
H0.bridge <- bridge_sampler(samples = samples_H0, data = data_H0,
log_posterior = log_posterior_H0, lb = lb_H0,
ub = ub_H0, silent = TRUE)
print(H0.bridge)

# compute log marginal likelihood via bridge sampling for H1
H1.bridge <- bridge_sampler(samples = samples_H1, data = data_H1,
log_posterior = log_posterior_H1, lb = lb_H1,
ub = ub_H1, silent = TRUE)
print(H1.bridge)

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

## ------------------------------------------------------------------------
# compute percentage errors
print(error_measures(H0.bridge)$percentage)
print(error_measures(H1.bridge)$percentage)
print(H0.error)
print(H1.error)

## ------------------------------------------------------------------------
# compute Bayes factor
Expand Down
33 changes: 22 additions & 11 deletions inst/doc/bridgesampling_example_jags.Rmd
View file
Expand Up @@ -37,7 +37,7 @@ y <- rnorm(n, theta, sqrt(sigma2))

Next, we specify the prior parameters $\mu_0$, $\tau^2_0$, $\alpha$, and $\beta$:

```{r}
```{r,eval=FALSE}
### set prior parameters ###
mu0 <- 0
tau20 <- 1
Expand All @@ -47,7 +47,7 @@ beta <- 1

## Fitting the Models
Now we can fit the null and the alternative model in `JAGS` (note that it is necessary to install `JAGS` for this). One usually requires a larger number of posterior sample for estimating the marginal likelihood than for simply estimating the model parameters. This is the reason for using a comparatively large number of samples (i.e., 50,000 post burn-in samples per chain) for this comparatively simple model.
```{r, message = FALSE, results='hide'}
```{r, eval=FALSE}
library(R2jags)
### functions to get posterior samples ###
Expand Down Expand Up @@ -118,7 +118,7 @@ When using MCMC software such as `JAGS` or `Stan`, specifying this function is r
The log of the densities on the right-hand side of these "`~`" symbols needs to be evaluated for the relevant quantities and then these log densities values are summed.

For example, in the null model, there are three "`~`" signs. Starting at the data-level, we need to evaluate the log of the normal density with mean $\theta_i$ and variance $\sigma^2$ for all $y_i$ and then sum the resulting log density values. Next, we move one step up in the model and evaluate the log of the group-level density for all $\theta_i$. Hence, we evaluate the log of the normal density for $\theta_i$ with mean $\mu = 0$ and variance $\tau^2$ (remember that `JAGS` parametrizes the normal distribution in terms of the precision `invTau2` = $1/\tau^2$; in contrast, `R` parametrizes it in terms of the standard deviation) and sum the resulting log density values. The result of this summation is added to the result of the previous summation for the data-level normal distribution. Finally, we need to evaluate the log of the prior density for `invTau2`. This means that we compute the log density of the gamma distribution with parameters $\alpha$ and $\beta$ for the sampled `invTau2` value and add the resulting log density value to the result of summing the data-level and group-level log densities. The unnormalized log posterior for the alternative model can be obtained in a similar fashion. The resulting functions look as follows:
```{r}
```{r,eval=FALSE}
### functions for evaluating the unnormalized posteriors on log scale ###
log_posterior_H0 <- function(samples.row, data) {
Expand Down Expand Up @@ -150,7 +150,7 @@ log_posterior_H1 <- function(samples.row, data) {

## Specifying the Parameter Bounds
The final step before computing the log marginal likelihoods is to specify the parameter bounds. In this example, for both models, all parameters can range from $-\infty$ to $\infty$ except the precision `invTau2` which has a lower bound of zero. These boundary vectors need to be named and the names need to match the order of the parameters.
```{r}
```{r,eval=FALSE}
# specify parameter bounds H0
cn <- colnames(samples_H0$BUGSoutput$sims.matrix)
cn <- cn[cn != "deviance"]
Expand All @@ -171,27 +171,38 @@ Note that currently, the lower and upper bound of a parameter cannot be a functi

## Computing the (Log) Marginal Likelihoods
Now we are ready to compute the log marginal likelihoods using the `bridge_sampler` function. We use `silent = TRUE` to suppress printing the number of iterations to the console:
```{r}
```{r, echo=FALSE}
load(system.file("extdata/", "vignette_example_jags.RData",
package = "bridgesampling"))
```

```{r,eval=FALSE}
# compute log marginal likelihood via bridge sampling for H0
H0.bridge <- bridge_sampler(samples = samples_H0, data = data_H0,
log_posterior = log_posterior_H0, lb = lb_H0,
ub = ub_H0, silent = TRUE)
print(H0.bridge)
# compute log marginal likelihood via bridge sampling for H1
H1.bridge <- bridge_sampler(samples = samples_H1, data = data_H1,
log_posterior = log_posterior_H1, lb = lb_H1,
ub = ub_H1, silent = TRUE)
```
We obtain:
```{r}
print(H0.bridge)
print(H1.bridge)
```

We can use the `error_measures` function to compute an approximate percentage error of the estimates:
```{r}
```{r,eval=FALSE}
# compute percentage errors
print(error_measures(H0.bridge)$percentage)
print(error_measures(H1.bridge)$percentage)
H0.error <- error_measures(H0.bridge)$percentage
H1.error <- error_measures(H1.bridge)$percentage
```
We obtain:
```{r}
print(H0.error)
print(H1.error)
```

## Bayesian Model Comparison
To compare the null model and the alternative model, we can compute the Bayes factor by using the `bf` function.
In our case, we compute $\text{BF}_{01}$, that is, the Bayes factor which quantifies how much more likely the data are under the null versus the alternative model:
Expand Down

0 comments on commit af8990f

Please sign in to comment.