CMCMC Example Workflows

What Is CMCMC?

CMCMC stands for Contemporaneous Markov chain Monte Carlo. It is designed for targets that are known up to a normalising constant,

\[ \pi(\theta) \propto f(\theta), \]

where the aim is to sample from \(\pi\) or estimate expectations such as \(E_\pi[h(\theta)]\). Standard MCMC can be slow when the target is highly correlated, high-dimensional, or difficult to explore. CMCMC addresses this by running many Metropolis chains in parallel and using the current population of chains to learn a useful proposal shape.

The main idea is to share information across chains at the same iteration. The population is split into two groups, \(A\) and \(B\). At a covariance refresh step, group \(A\) uses a covariance estimated from the current states in group \(B\), and group \(B\) uses a covariance estimated from group \(A\):

\[ \widehat{C}_{t,A} = \operatorname{Cov}\{\theta^{(k)}_t : k \in A\}, \qquad \widehat{C}_{t,B} = \operatorname{Cov}\{\theta^{(k)}_t : k \in B\}, \]

For example, a chain in group \(A\) proposes

\[ \theta^\star = \theta^{(k)}_t + z, \qquad z \sim N_D(0, s_D\widehat{C}_{t,B}), \qquad s_D = \frac{2.38^2}{D}, \]

and then accepts the proposal with the usual Metropolis probability \(\min\{1, f(\theta^\star)/f(\theta^{(k)}_t)\}\). The update is still a Metropolis step for the same target; the difference is how the proposal covariance is learned. Conditional on the covariance borrowed from the other group, the proposal is symmetric, so the usual Metropolis ratio is unchanged. The same scaling \(s_D\) is used by the CPU and CUDA backends. Since each group contains \(K/2\) particles, the package requires \(K/2 > D\). This avoids the immediate rank-deficient case where an empirical covariance matrix cannot be full rank in \(D\) dimensions.

This package was developed to make that idea practical on modern parallel hardware. Instead of relying on a long history from one chain, CMCMC uses the current population of \(K\) chains. This fits GPU execution well, because many target densities can be evaluated at the same time, and the sampler only needs to keep the current \(K \times D\) state array in memory rather than a full trajectory history.

This vignette gives two examples of using the package: a multivariate normal target and a real-data hierarchical model with a user-supplied target density. The examples use the regular step mode so that iteration traces can be saved and plotted. A separate section below shows how to switch to CUDA batch mode when only the final population is needed. GLM examples using glm_cmcmc() are covered in the separate GLM CMCMC vignette.

The examples use backend = "auto". With this setting, CMCMC uses CUDA when nvcc is available and otherwise falls back to the CPU backend.

When the vignette is rendered, the sampler sections produce summary tables and plots. On machines without CUDA, the same examples run on the CPU with smaller settings so that they remain practical.

Quick Start

The smallest useful call only needs an initial population, a target kernel, the number of iterations, the covariance-refresh interval, and the data vectors required by the kernel. The example below uses the built-in multivariate normal kernel and the CPU backend, so it can run on machines without CUDA.

library(CMCMC)

D <- 2L
K <- 100L
set.seed(1)

mu <- c(0.25, -0.5)
Sigma <- matrix(c(1, 0.4, 0.4, 1), nrow = D)
X <- c(mu, as.numeric(solve(Sigma)))
Xi <- integer(0)
init <- matrix(rnorm(K * D), nrow = K, ncol = D)

samples <- CMCMC::cmcmc(
  init = init,
  GPUkernel = "MVNorm",
  it = 50,
  covit = 5,
  k = K,
  X = X,
  Xi = Xi,
  saved_iterations = 1,
  seed = 123,
  backend = "cpu"
)

dim(samples)
#> [1] 100   4
colMeans(samples[, 3:(2 + D), drop = FALSE])
#> [1]  0.05521669 -0.73397789

Backend Selection

The package also includes inca(), which implements interchain adaptive Markov chain Monte Carlo (INCA) in its global covariance form. This follows the parallel-chain adaptation idea of Craiu, Rosenthal and Yang (2009), where proposal covariance information is pooled across chains after burn-in.

Both cmcmc() and inca() accept backend = "auto", backend = "cuda", and backend = "cpu".

The default auto mode is the most convenient option. It uses CUDA when a CUDA compiler is available and uses the CPU backend otherwise.

Use backend = "cuda" when CUDA is required and the call should fail if CUDA is not available. Use backend = "cpu" when the run should stay on the CPU. The CPU backend is always built. If OpenMP is available at installation time it can use multiple CPU threads through nthreads; otherwise the same CPU backend runs in serial mode. If CUDA is not available at installation time, the package still installs with a CUDA stub backend. In that case backend = "cpu" works, while backend = "cuda" reports that CUDA is unavailable. For the CPU backend, random numbers are generated serially with R’s RNG before the OpenMP particle-update loop. This keeps R’s global RNG state out of the parallel region while still parallelising the chain updates.

CMCMC tries to detect CUDA paths when the package is loaded. If it does not find your CUDA installation, set the CUDA binary and library directories manually before using backend = "cuda" or a user-provided CUDA kernel. The cmcmcopt() settings are session-local. If you want these paths to be used automatically in future R sessions, set CMCMC_CUDA_BIN and CMCMC_CUDA_LIB64 before loading CMCMC, for example in ~/.Renviron.

library(CMCMC)

CMCMC::cmcmcopt(
  PATH = "/usr/local/cuda/bin",
  DYLD_LIBRARY_PATH = "/usr/local/cuda/lib64"
)

On Windows, DYLD_LIBRARY_PATH is used by CMCMC as the CUDA library directory and is also prepended to PATH when compiling runtime CUDA plugins.

CUDA target kernels are compiled at runtime into a user-writable cache. By default this is under file.path(tools::R_user_dir("CMCMC", "cache"), "kernels"). If that location is not suitable, set CMCMC_KERNEL_CACHE to a writable kernel-cache root before running CUDA examples. User CUDA kernels installed with load_target_kernel() are stored in the same kernel cache, so the package does not need write access to its installed library directory. Launch-tuning results are cached separately under tools::R_user_dir("CMCMC", "cache"). To remove runtime kernel files, for example before uninstalling the package, call cmcmc_clear_kernel_cache().

Step Mode, Batch Mode, and Launch Tuning

The default CUDA sampler uses step mode. In step mode the package computes a proposal covariance every covit iterations, generates proposals with CUDA libraries, and launches the Metropolis kernel once per iteration until the next covariance refresh. This mode can save multiple iterations, so it is useful for trace plots and diagnostics.

With batch = TRUE, the CUDA backend pre-generates the proposal and uniform random numbers for one covariance-refresh block, then runs the covit Metropolis updates inside one kernel call. This reduces repeated kernel launches and intermediate global-memory writes. Batch mode returns only the final population, so saved_iterations is set to 1. Batch mode is only used by the CUDA backend; the CPU backend uses the regular sampler.

samples <- CMCMC::cmcmc(
  init = init,
  GPUkernel = "MVNorm",
  it = 1000,
  covit = 20,
  k = K,
  X = X,
  Xi = Xi,
  saved_iterations = 1,
  seed = 1,
  backend = "cuda",
  batch = TRUE
)

The CUDA backend chooses a launch configuration automatically. If no tuned configuration is available, it asks a CUDA occupancy helper for a default block size. For a specific target, dimension, particle count, covit, and batch setting, cmcmc_autotune() can run a short low-level timing experiment and store the fastest launch configuration in the user cache. It returns a sampler function that reuses the tuned launch.

tuned_sampler <- CMCMC::cmcmc_autotune(
  init = init,
  GPUkernel = "MVNorm",
  covit = 20,
  k = K,
  X = X,
  Xi = Xi,
  batch = TRUE,
  repeats = 3
)

attr(tuned_sampler, "launch")
samples <- tuned_sampler(it = 1000, seed = 1)
CMCMC::cmcmc_autotune_list()

Common Output Layout

Both cmcmc() and inca() return a matrix. The first column is the iteration, the second column is the particle index, and the remaining columns are the sampled state vector. Here a particle means one member of the parallel chain population.

theta <- samples[, 3:ncol(samples), drop = FALSE]
posterior_mean <- colMeans(theta)

When saved_iterations = 0, all iterations are saved. This is useful for diagnostic plots, but larger runs should usually save only the final part of the run to keep output size manageable.

extract_iter_mean_trace <- function(samples, d) {
  iter <- samples[, 1]
  vals <- samples[, 3:(2 + d), drop = FALSE]
  sums <- rowsum(vals, group = iter, reorder = TRUE)
  counts <- as.numeric(table(iter)[rownames(sums)])
  means <- sweep(sums, 1, counts, "/")
  list(iter = as.integer(rownames(sums)), mean = means)
}

Multivariate Normal

The MVNorm kernel evaluates a multivariate normal target density.

The numeric data vector is X = c(mu, SigmaInv), where SigmaInv stores \(\Sigma^{-1}\) in R’s column-major order. Xi is empty.

In this example there is no observed data model. We specify the target distribution directly. Let \(\theta\) be a \(D\)-dimensional real-valued vector. Then

\[ \theta \sim N_D(\mu, \Sigma), \]

with target density

\[ \pi(\theta) \propto |\Sigma|^{-1/2} \exp\!\left\{-\frac{1}{2}(\theta-\mu)^\top\Sigma^{-1}(\theta-\mu)\right\}. \]

library(CMCMC)

sampler <- "cmcmc"
D <- 10L
K <- if (has_nvcc) 1000L else 200L
it <- if (has_nvcc) 1000L else 300L
covit <- 10L
rho <- 0.9
seed <- 1L

set.seed(seed)

mu <- rnorm(D, mean = 0, sd = 1)
Sigma <- matrix(rho, nrow = D, ncol = D)
diag(Sigma) <- 1

X <- c(mu, as.numeric(solve(Sigma)))
Xi <- integer(0)
init <- matrix(runif(K * D, -5, 6), nrow = K, ncol = D)

sampler_fun <- if (sampler == "inca") CMCMC::inca else CMCMC::cmcmc
samples <- sampler_fun(
  init = init,
  GPUkernel = "MVNorm",
  it = it,
  covit = covit,
  k = K,
  X = X,
  Xi = Xi,
  saved_iterations = 0,
  seed = seed,
  backend = "auto"
)

theta <- samples[, 3:(2 + D), drop = FALSE]
last_iter <- max(samples[, 1])
tail_iters <- min(200L, last_iter)
theta_tail <- samples[samples[, 1] > (last_iter - tail_iters),
                      3:(2 + D), drop = FALSE]

mvnorm_summary <- data.frame(
  dimension = paste0("theta[", seq_len(D), "]"),
  posterior_mean = colMeans(theta_tail),
  target_mean = mu,
  row.names = NULL
)

mvnorm_summary
#>    dimension posterior_mean target_mean
#> 1   theta[1]     -0.6103807  -0.6264538
#> 2   theta[2]      0.2068573   0.1836433
#> 3   theta[3]     -0.7924175  -0.8356286
#> 4   theta[4]      1.6095096   1.5952808
#> 5   theta[5]      0.3316905   0.3295078
#> 6   theta[6]     -0.8226705  -0.8204684
#> 7   theta[7]      0.5214049   0.4874291
#> 8   theta[8]      0.7585787   0.7383247
#> 9   theta[9]      0.6022115   0.5757814
#> 10 theta[10]     -0.2865521  -0.3053884

The summary table compares posterior means from the saved tail of the run with the known target mean. This gives a direct check that the sampler is centered in the right region.

The first plot shows the population mean trace for the first five dimensions. The line is the mean across particles at each iteration, and the dashed line is the true target mean.

The second plot shows marginal distributions for the same dimensions, using the last saved iterations pooled across particles. The dashed curve is the corresponding target marginal density.

trace <- extract_iter_mean_trace(samples, D)
dims_to_plot <- seq_len(min(5L, D))

op <- par(no.readonly = TRUE)
par(mfrow = c(length(dims_to_plot), 1))

for (j in dims_to_plot) {
  plot(
    trace$iter, trace$mean[, j],
    type = "l",
    xlab = "Iteration",
    ylab = paste0("mean(theta[", j, "])")
  )
  abline(h = mu[j], lty = 2)
}


par(op)

op <- par(no.readonly = TRUE)
par(mfrow = c(2, 3))

for (j in dims_to_plot) {
  hist(
    theta_tail[, j],
    breaks = 35,
    probability = TRUE,
    main = paste0("theta[", j, "]"),
    xlab = "Sample value"
  )
  lines(density(theta_tail[, j]))
  curve(
    dnorm(x, mean = mu[j], sd = sqrt(Sigma[j, j])),
    add = TRUE,
    lty = 2
  )
}

plot.new()
legend(
  "center",
  legend = c("sample density", "target marginal"),
  lty = c(1, 2)
)


par(op)

Hierarchical Beta-Binomial Rat Tumor Rates

When the target density is not one of the built-in kernels, users can provide their own kernel. There are two routes:

load_target_kernel() installs a user CUDA .cu kernel into CMCMC’s user cache for the CUDA backend.
load_target_kernel_cpu() compiles a user C .c kernel for the CPU backend. The compiled CPU shared library is created under tempdir(), so re-run load_target_kernel_cpu() in a new R session before using that CPU kernel again.

The example below uses the same hierarchical beta-binomial model in both routes.

If CUDA is available, the vignette uses the CUDA kernel. Otherwise, it compiles the plain C kernel and runs on the CPU. The sampler call still uses backend = "auto".

User Kernel Interface

The user-provided kernel is only responsible for evaluating the target log density. CMCMC still handles the population of chains, proposal generation, Metropolis accept/reject steps, covariance refreshes, and sample storage. This keeps the user kernel small: it receives one particle state theta, the numeric data vector X, and the integer data vector Xi, and it returns the log target density up to an additive constant.

The user decides the layout of X and Xi. The sampler only copies these vectors to the target-density function. A useful convention is to put continuous quantities in X and integer quantities in Xi. For example, X can contain a design matrix, observed numeric responses, or prior hyperparameters, while Xi can contain sample sizes, numbers of groups, and count data. The same layout must then be used inside logdens().

For the CUDA backend, the kernel file must define a device function with the signature

__device__ float logdens(const float *theta, const float *Xf, const int *Xi)

The package compiles this function into the CUDA Metropolis kernel. The symbols DIM and K are supplied at compile time by CMCMC, so the CUDA log-density can check that the state dimension matches the model. The vector theta is the current particle state, Xf contains numeric data copied to the GPU, and Xi contains integer data copied to the GPU.

For the CPU backend, the kernel file must define a C function with the signature

double logdens(const double *theta, int dim,
               const double *X, int lenX,
               const int *Xi, int lenXi)

The CPU version receives the same model information, but it also receives the dimension and data-vector lengths directly. In both versions, impossible parameter values should return -INFINITY; this makes the Metropolis step reject the proposal.

Model

The data are rat tumor counts from Tarone’s historical control experiments. They are reproduced in Bayesian Data Analysis and distributed with the book’s example data.

For experiment \(j=1,\ldots,J\), let \(y_j\) be the number of rats with tumors out of \(n_j\) rats. The model is

\[ y_j \mid \omega_j \sim \operatorname{Binomial}(n_j, \omega_j), \qquad \omega_j \mid \alpha,\beta \sim \operatorname{Beta}(\alpha,\beta). \]

The hyperprior is

\[ \pi(\alpha,\beta) \propto (\alpha+\beta)^{-5/2}. \]

Ignoring binomial coefficients that are constant in the parameters, the posterior is

\[ \pi(\alpha,\beta,\omega \mid y,n) \propto (\alpha+\beta)^{-5/2} \prod_{j=1}^{J} \left[ \omega_j^{y_j}(1-\omega_j)^{n_j-y_j} \frac{\omega_j^{\alpha-1}(1-\omega_j)^{\beta-1}} {B(\alpha,\beta)} \right]. \]

The user kernel below works with the unconstrained state (log(alpha), log(beta), z_1, ..., z_J), where z_j = logit(omega_j).

The integer data vector is Xi = c(J, n_1, ..., n_J, y_1, ..., y_J), and X is empty.

The kernel data layout is:

Object	Contents	Notes
`theta[1]`	`log(alpha)`	unconstrained hyperparameter
`theta[2]`	`log(beta)`	unconstrained hyperparameter
`theta[3:(2 + J)]`	`logit(omega_1), ..., logit(omega_J)`	unconstrained experiment rates
`X`	empty numeric vector	no continuous data are needed
`Xi`	`J`, then all `n_j`, then all `y_j`	integer data for the kernel

Because the sampler uses this unconstrained state, the kernel includes the Jacobian terms log(alpha) + log(beta) and log(omega_j) + log(1 - omega_j). With these terms included, the kernel targets the constrained posterior written above.

The CUDA and CPU kernels below therefore have the same structure. They read J, n_j, and y_j from Xi; transform the unconstrained parameters to alpha, beta, and omega_j; add the hyperprior, beta-binomial likelihood, beta-density contribution, and Jacobian adjustment; and return the resulting log density. The CUDA version uses float, CUDA device functions, and single precision math functions such as expf() and lgammaf(). The CPU version uses double and the corresponding C math functions.

This is not the normal treatment-effect rat example in Robert (2007, Example 10.2.3, Chapter 10). Robert’s example models control, intoxication, placebo, and drug effects with normal hierarchical levels.

The example here models binomial tumor rates. It is included because it is a compact real-data example for a user-supplied target density.

CUDA and CPU Kernels

rat_cuda_kernel_code <- '
#include <math.h>

static __device__ __forceinline__ float log_sigmoid(float x) {
  return -log1pf(expf(-x));
}

static __device__ __forceinline__ float log1m_sigmoid(float x) {
  return -log1pf(expf(x));
}

__device__ float logdens(const float *theta, const float *Xf, const int *Xi) {
  (void)Xf;
  if (!Xi) return -INFINITY;

  const int J = Xi[0];
  if (J <= 0 || DIM != (2 + J)) return -INFINITY;

  const float log_alpha = theta[0];
  const float log_beta = theta[1];
  const float alpha = expf(log_alpha);
  const float beta = expf(log_beta);
  if (!isfinite(alpha) || !isfinite(beta) || alpha <= 0.0f || beta <= 0.0f) {
    return -INFINITY;
  }

  const float s = alpha + beta;
  if (!isfinite(s) || s <= 0.0f) return -INFINITY;

  float lp = -2.5f * logf(s);
  lp += log_alpha + log_beta;

  const float logB_ab = lgammaf(alpha) + lgammaf(beta) - lgammaf(s);

  for (int j = 0; j < J; j++) {
    const int nj = Xi[1 + j];
    const int yj = Xi[1 + J + j];
    if (nj < 0 || yj < 0 || yj > nj) return -INFINITY;

    const float z = theta[2 + j];
    const float lomega = log_sigmoid(z);
    const float l1momega = log1m_sigmoid(z);

    lp += ((float)yj) * lomega + ((float)(nj - yj)) * l1momega;
    lp += (alpha - 1.0f) * lomega + (beta - 1.0f) * l1momega - logB_ab;
    lp += lomega + l1momega;
  }

  return lp;
}
'

rat_cpu_kernel_code <- '
#include <math.h>

#ifdef _WIN32
__declspec(dllexport)
#endif
double logdens(const double *theta, int dim,
               const double *X, int lenX,
               const int *Xi, int lenXi) {
  (void)X;
  (void)lenX;
  if (!Xi || lenXi < 1) return -INFINITY;

  const int J = Xi[0];
  if (J <= 0 || dim != (2 + J) || lenXi < 1 + 2 * J) return -INFINITY;

  const double log_alpha = theta[0];
  const double log_beta = theta[1];
  const double alpha = exp(log_alpha);
  const double beta = exp(log_beta);
  if (!isfinite(alpha) || !isfinite(beta) || alpha <= 0.0 || beta <= 0.0) {
    return -INFINITY;
  }

  const double s = alpha + beta;
  if (!isfinite(s) || s <= 0.0) return -INFINITY;

  double lp = -2.5 * log(s) + log_alpha + log_beta;
  const double logB_ab = lgamma(alpha) + lgamma(beta) - lgamma(s);

  for (int j = 0; j < J; j++) {
    const int nj = Xi[1 + j];
    const int yj = Xi[1 + J + j];
    if (nj < 0 || yj < 0 || yj > nj) return -INFINITY;

    const double z = theta[2 + j];
    const double lomega = (z >= 0.0) ? -log1p(exp(-z)) : z - log1p(exp(z));
    const double l1momega = (z >= 0.0) ? -z - log1p(exp(-z)) : -log1p(exp(z));

    lp += ((double)yj) * lomega + ((double)(nj - yj)) * l1momega;
    lp += (alpha - 1.0) * lomega + (beta - 1.0) * l1momega - logB_ab;
    lp += lomega + l1momega;
  }

  return lp;
}
'

if (has_nvcc) {
  rat_kernel_path <- tempfile(fileext = ".cu")
  writeLines(rat_cuda_kernel_code, rat_kernel_path)
  rat_kernel <- CMCMC::load_target_kernel(
    rat_kernel_path,
    name = "RatTumorBetaBinom.cu"
  )
} else {
  rat_kernel_path <- tempfile(fileext = ".c")
  writeLines(rat_cpu_kernel_code, rat_kernel_path)
  rat_kernel <- CMCMC::load_target_kernel_cpu(
    rat_kernel_path,
    name = "RatTumorBetaBinom"
  )
}
#> Compiling CPU plugin: kernel=RatTumorBetaBinom
#>   output: /tmp/RtmpjA8R9D/CMCMC_cpu_kernels/RatTumorBetaBinom/RatTumorBetaBinom.so

Running the Sampler

library(CMCMC)

rat_tumor <- as.data.frame(matrix(c(
  0, 20, 0, 20, 0, 20, 0, 20, 0, 20, 0, 20, 0, 20,
  0, 19, 0, 19, 0, 19, 0, 19, 0, 18, 0, 18, 0, 17,
  1, 20, 1, 20, 1, 20, 1, 20, 1, 19, 1, 19, 1, 18,
  1, 18, 3, 27, 2, 25, 2, 24, 2, 23, 2, 20, 2, 20,
  2, 20, 2, 20, 2, 20, 2, 20, 1, 10, 5, 49, 2, 19,
  5, 46, 2, 17, 7, 49, 7, 47, 3, 20, 3, 20, 2, 13,
  9, 48, 10, 50, 4, 20, 4, 20, 4, 20, 4, 20, 4, 20,
  4, 20, 4, 20, 10, 48, 4, 19, 4, 19, 4, 19, 5, 22,
  11, 46, 12, 49, 5, 20, 5, 20, 6, 23, 5, 19, 6, 22,
  6, 20, 6, 20, 6, 20, 16, 52, 15, 46, 15, 47, 9, 24,
  4, 14
), ncol = 2, byrow = TRUE))
names(rat_tumor) <- c("y", "n")

y <- as.integer(rat_tumor$y)
n <- as.integer(rat_tumor$n)
J <- length(y)
D <- 2L + J
K <- if (has_nvcc) 512L else 256L
it <- if (has_nvcc) 2000L else 1200L
covit <- 20L
saved_iterations <- if (has_nvcc) 500L else 300L
seed <- 1L

Xi <- as.integer(c(J, n, y))
X <- numeric(0)

set.seed(seed)
p0 <- (y + 0.5) / (n + 1.0)
p0 <- pmin(pmax(p0, 1e-6), 1 - 1e-6)
init_center <- c(log(1.4), log(8.6), qlogis(p0))
init <- matrix(rnorm(K * D, mean = init_center, sd = 0.15),
               nrow = K, ncol = D, byrow = TRUE)

samples <- CMCMC::cmcmc(
  init = init,
  GPUkernel = rat_kernel,
  it = it,
  covit = covit,
  k = K,
  X = X,
  Xi = Xi,
  saved_iterations = saved_iterations,
  seed = seed,
  backend = "auto"
)

theta <- samples[, 3:(2 + D), drop = FALSE]
alpha <- exp(theta[, 1])
beta <- exp(theta[, 2])
mu <- alpha / (alpha + beta)
size <- alpha + beta
omega <- plogis(theta[, 3:ncol(theta), drop = FALSE])

rat_hyper_summary <- data.frame(
  quantity = c(
    "alpha",
    "beta",
    "population mean",
    "alpha + beta"
  ),
  posterior_mean = c(
    mean(alpha),
    mean(beta),
    mean(mu),
    mean(size)
  ),
  posterior_median = c(
    median(alpha),
    median(beta),
    median(mu),
    median(size)
  ),
  row.names = NULL
)

rat_hyper_summary
#>          quantity posterior_mean posterior_median
#> 1           alpha      2.3729052         2.238163
#> 2            beta     14.0778646        13.340850
#> 3 population mean      0.1446606         0.144110
#> 4    alpha + beta     16.4507699        15.556339

log_sum_exp <- function(x) {
  m <- max(x)
  m + log(sum(exp(x - m)))
}

rat_marginal_grid <- function(y, n, n_eta = 260L, n_log_size = 260L) {
  eta_grid <- seq(-3.5, -0.4, length.out = n_eta)
  log_size_grid <- seq(log(1), log(80), length.out = n_log_size)
  grid <- expand.grid(
    eta = eta_grid,
    log_size = log_size_grid
  )

  grid$mu <- plogis(grid$eta)
  grid$size <- exp(grid$log_size)
  grid$alpha <- grid$size * grid$mu
  grid$beta <- grid$size * (1 - grid$mu)

  log_post <- -2.5 * log(grid$size)
  for (j in seq_along(y)) {
    log_post <- log_post +
      lbeta(grid$alpha + y[j], grid$beta + n[j] - y[j]) -
      lbeta(grid$alpha, grid$beta)
  }

  # The grid is uniform in (log(alpha / beta), log(alpha + beta)).
  # Add the Jacobian for the transformation back to (alpha, beta).
  log_weight <- log_post + log(grid$alpha) + log(grid$beta)
  grid$weight <- exp(log_weight - log_sum_exp(log_weight))
  grid
}

rat_grid <- rat_marginal_grid(y, n)

grid_reference <- c(
  sum(rat_grid$weight * rat_grid$alpha),
  sum(rat_grid$weight * rat_grid$beta),
  sum(rat_grid$weight * rat_grid$mu),
  sum(rat_grid$weight * rat_grid$size)
)

bda_moment_alpha <- 1.4
bda_moment_beta <- 8.6
bda_moment_size <- bda_moment_alpha + bda_moment_beta

bda_moment_reference <- c(
  bda_moment_alpha,
  bda_moment_beta,
  bda_moment_alpha / bda_moment_size,
  bda_moment_size
)

bda_moment_prior_summary <- data.frame(
  quantity = rat_hyper_summary$quantity,
  bda_moment_estimate = bda_moment_reference,
  row.names = NULL
)

bda_moment_prior_summary
#>          quantity bda_moment_estimate
#> 1           alpha                1.40
#> 2            beta                8.60
#> 3 population mean                0.14
#> 4    alpha + beta               10.00

rat_reference_summary <- data.frame(
  quantity = rat_hyper_summary$quantity,
  cmcmc_posterior_mean = rat_hyper_summary$posterior_mean,
  grid_posterior_mean = grid_reference,
  absolute_difference = abs(rat_hyper_summary$posterior_mean - grid_reference),
  row.names = NULL
)

rat_reference_summary
#>          quantity cmcmc_posterior_mean grid_posterior_mean absolute_difference
#> 1           alpha            2.3729052            2.401469        0.0285635539
#> 2            beta           14.0778646           14.313147        0.2352819431
#> 3 population mean            0.1446606            0.144297        0.0003635232
#> 4    alpha + beta           16.4507699           16.714615        0.2638454970

current_experiment <- J
bda_current_alpha <- bda_moment_alpha + y[current_experiment]
bda_current_beta <- bda_moment_beta + n[current_experiment] - y[current_experiment]
bda_current_rate_mean <- bda_current_alpha / (bda_current_alpha + bda_current_beta)

grid_current_rate_mean <- sum(
  rat_grid$weight *
    ((rat_grid$alpha + y[current_experiment]) /
       (rat_grid$alpha + rat_grid$beta + n[current_experiment]))
)

current_rate_comparison <- data.frame(
  method = c(
    "raw current experiment",
    "BDA moment-prior posterior",
    "grid marginal posterior",
    "CMCMC hierarchical posterior"
  ),
  estimate = c(
    y[current_experiment] / n[current_experiment],
    bda_current_rate_mean,
    grid_current_rate_mean,
    mean(omega[, current_experiment])
  ),
  row.names = NULL
)

current_rate_comparison
#>                         method  estimate
#> 1       raw current experiment 0.2857143
#> 2   BDA moment-prior posterior 0.2250000
#> 3      grid marginal posterior 0.2108613
#> 4 CMCMC hierarchical posterior 0.2058705

Posterior Checks

Gelman et al. (2013) use the historical rat experiments to motivate the moment estimate \((\alpha,\beta)=(1.4,8.6)\) for the beta population distribution.

That value is useful context, but it is not the posterior target sampled here. If it is used as a fixed prior for the current \(4/14\) experiment, the conjugate posterior is \(\operatorname{Beta}(5.4,18.6)\), with posterior mean \(0.225\).

For the posterior check, we use a deterministic grid approximation to \(\pi(\alpha,\beta\mid y,n)\). The grid calculation integrates out the experiment-specific rates \(\omega_j\).

That reference posterior is proportional to

\[ \pi(\alpha,\beta\mid y,n) \propto (\alpha+\beta)^{-5/2} \prod_{j=1}^{J} \frac{B(\alpha+y_j,\beta+n_j-y_j)}{B(\alpha,\beta)}. \]

The rat_reference_summary table compares CMCMC posterior means with this grid reference. Small differences are expected because CMCMC uses Monte Carlo sampling and the grid is only a numerical approximation.

For the current experiment, where \(y=4\) out of \(n=14\) rats developed tumors, the current_rate_comparison table compares the raw rate, the grid posterior mean, and the CMCMC posterior mean for the same experiment.

The population mean tumor rate is \(\alpha/(\alpha+\beta)\).

The quantity \(\alpha+\beta\) controls how tightly the experiment-specific rates are pulled toward that population mean. Larger values mean stronger pooling. Smaller values allow more variation between experiments.

The experiment-specific posterior tumor rates shrink the raw rates \(y_j/n_j\) toward the population distribution learned from all experiments.

The next table reports the first ten experiment-specific rates.

Conditional on \((\alpha,\beta)\), the posterior mean for one experiment-specific rate is

\[ E(\omega_j \mid y_j,n_j,\alpha,\beta) = \frac{\alpha+y_j}{\alpha+\beta+n_j}. \]

This is why the posterior means are shrunk toward the population mean. Experiments with \(y_j=0\) are not estimated as exactly zero, and experiments with high raw rates are pulled downward.

The uncertainty intervals are wider when the experiment sample size is small or the raw rate is far from the overall pattern.

omega_summary <- data.frame(
  experiment = seq_len(J),
  tumors = y,
  rats = n,
  raw_rate = y / n,
  posterior_mean = colMeans(omega),
  q05 = apply(omega, 2, quantile, probs = 0.05),
  q95 = apply(omega, 2, quantile, probs = 0.95),
  row.names = NULL
)

head(omega_summary, 10)
#>    experiment tumors rats raw_rate posterior_mean        q05       q95
#> 1           1      0   20        0     0.06344772 0.01289059 0.1404055
#> 2           2      0   20        0     0.06119992 0.01014063 0.1400613
#> 3           3      0   20        0     0.06204706 0.01052461 0.1375315
#> 4           4      0   20        0     0.06375070 0.01245682 0.1421737
#> 5           5      0   20        0     0.06238221 0.01150523 0.1376593
#> 6           6      0   20        0     0.06466333 0.01157556 0.1503766
#> 7           7      0   20        0     0.06205197 0.01085246 0.1448190
#> 8           8      0   19        0     0.06616161 0.01043631 0.1484563
#> 9           9      0   19        0     0.06628421 0.01368820 0.1473706
#> 10         10      0   19        0     0.06712421 0.01008265 0.1520940

The next figure shows simple diagnostics for the hyperparameters. The first three panels track particle means for \(\alpha\), \(\beta\), and the population mean tumor rate. The fourth panel shows the posterior density of \(\alpha/(\alpha+\beta)\).

rat_trace <- extract_iter_mean_trace(
  cbind(samples[, 1:2, drop = FALSE], alpha = alpha, beta = beta, mu = mu, size = size),
  4L
)

op <- par(no.readonly = TRUE)
par(mfrow = c(2, 2))
plot(rat_trace$iter, rat_trace$mean[, 1], type = "l",
     xlab = "Iteration", ylab = "mean(alpha)", main = "alpha")
plot(rat_trace$iter, rat_trace$mean[, 2], type = "l",
     xlab = "Iteration", ylab = "mean(beta)", main = "beta")
plot(rat_trace$iter, rat_trace$mean[, 3], type = "l",
     xlab = "Iteration", ylab = "mean population rate", main = "alpha/(alpha+beta)")
hist(mu, breaks = 35, probability = TRUE,
     main = "Population rate posterior", xlab = "alpha/(alpha+beta)")
lines(density(mu))

par(op)

The final rat plot compares each raw rate with its posterior mean and 90% posterior interval. The experiments are ordered by raw rate, which makes the shrinkage pattern easier to see.

ord <- order(omega_summary$raw_rate)
plot(
  omega_summary$raw_rate[ord],
  pch = 19,
  col = "gray55",
  ylim = range(c(omega_summary$raw_rate, omega_summary$q05, omega_summary$q95)),
  xlab = "Experiment ordered by raw rate",
  ylab = "Tumor rate",
  main = "Rat tumor rates: raw estimates and posterior shrinkage"
)
segments(
  x0 = seq_along(ord),
  y0 = omega_summary$q05[ord],
  x1 = seq_along(ord),
  y1 = omega_summary$q95[ord],
  col = "gray75"
)
points(omega_summary$posterior_mean[ord], pch = 19)
legend(
  "topleft",
  legend = c("raw y/n", "posterior mean", "90% interval"),
  pch = c(19, 19, NA),
  lty = c(NA, NA, 1),
  col = c("gray55", "black", "gray75"),
  bty = "n"
)

References

Craiu, R. V., Rosenthal, J. S. and Yang, C. (2009). Learn From Thy Neighbor: Parallel-Chain and Regional Adaptive MCMC. Journal of the American Statistical Association, 104(488), 1454-1466. doi:10.1198/jasa.2009.tm08393.

Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A. and Rubin, D. B. (2013). Bayesian Data Analysis. 3rd edition. Chapman and Hall/CRC.

Gelman, A., Carlin, J. B., Stern, H. S. and Rubin, D. B. (2003). Bayesian Data Analysis. 2nd edition. Chapman and Hall/CRC. Rat tumor example data: https://sites.stat.columbia.edu/gelman/book/data/rats.asc.

Robert, C. P. (2007). The Bayesian Choice: From Decision-Theoretic Foundations to Computational Implementation. 2nd edition. Springer.

Tarone, R. E. (1982). The use of historical control information in testing for a trend in proportions. Biometrics, 38, 215-220.