---
title: "CMCMC Example Workflows"
output:
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{CMCMC Example Workflows}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE, purl=FALSE}
nvcc_name <- if (.Platform$OS.type == "windows") "nvcc.exe" else "nvcc"
has_nvcc <- nzchar(Sys.which(nvcc_name))
run_cmcmc_examples <- TRUE

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = FALSE,
  purl = FALSE
)
```

## 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.

```{r quick-start, eval=exists("run_cmcmc_examples") && run_cmcmc_examples, purl=FALSE}
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)
colMeans(samples[, 3:(2 + D), drop = FALSE])
```

## 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`.

```{r cuda-setup, purl=FALSE}
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.

```{r batch-mode-example, eval=FALSE, purl=FALSE}
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.

```{r tuning-example, eval=FALSE, purl=FALSE}
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.

```{r output-layout, purl=FALSE}
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.

```{r trace-helper, eval=exists("run_cmcmc_examples") && run_cmcmc_examples, purl=FALSE}
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\}.
\]

```{r mvnorm, eval=exists("run_cmcmc_examples") && run_cmcmc_examples, purl=FALSE}
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
```

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.

```{r mvnorm-trace-plot, eval=exists("run_cmcmc_examples") && run_cmcmc_examples, fig.width=7, fig.height=8, purl=FALSE}
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)
```

```{r mvnorm-marginal-distribution, eval=exists("run_cmcmc_examples") && run_cmcmc_examples, fig.width=7, fig.height=5, purl=FALSE}
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

```c
__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

```c
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](https://sites.stat.columbia.edu/gelman/book/data/rats.asc).

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

```{r rat-tumor-user-kernel, eval=exists("run_cmcmc_examples") && run_cmcmc_examples, results='hide', purl=FALSE}
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"
  )
}
```

### Running the Sampler

```{r rat-tumor-data-and-sampling, eval=exists("run_cmcmc_examples") && run_cmcmc_examples, purl=FALSE}
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

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

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

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
```

### 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.

```{r rat-tumor-rate-summary, eval=exists("run_cmcmc_examples") && run_cmcmc_examples, purl=FALSE}
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)
```

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)\).

```{r rat-tumor-hyperparameter-plot, eval=exists("run_cmcmc_examples") && run_cmcmc_examples, fig.width=7, fig.height=5, purl=FALSE}
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.

```{r rat-tumor-rate-plot, eval=exists("run_cmcmc_examples") && run_cmcmc_examples, fig.width=7, fig.height=5, purl=FALSE}
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.