/computational-inference

name: computational-inference
description: Computational methods for statistical inference and optimization

Computational Inference

**Advanced computational methods for statistical inference in complex models**

Use this skill when working on: MCMC algorithms, importance sampling, Bayesian inference, parallel computing for statistics, GPU acceleration, or computationally intensive inference procedures.

---

Monte Carlo Methods

Fundamental Principle

Monte Carlo methods approximate expectations via random sampling:

$$E[g(X)] \approx \frac{1}{N} \sum_{i=1}^{N} g(X_i), \quad X_i \sim P$$

**Monte Carlo Standard Error**: $$\text{MCSE} = \frac{\hat{\sigma}}{\sqrt{N}}$$

Variance Reduction Techniques

| Technique | Idea | Variance Reduction | |-----------|------|-------------------| | Antithetic variates | Use negatively correlated pairs | Up to 50% | | Control variates | Subtract known expectation | Depends on correlation | | Importance sampling | Sample from better distribution | Can be dramatic | | Stratified sampling | Sample from strata separately | Reduces variance |

R Implementation

#' Monte Carlo Integration with Variance Reduction
#'
#' @param g Function to integrate
#' @param sampler Function that generates samples
#' @param n Number of samples
#' @param method Variance reduction method
#' @return Estimate with standard error
monte_carlo_integrate <- function(g, sampler, n = 10000,
                                   method = c("naive", "antithetic", "control")) {
  method <- match.arg(method)

  if (method == "naive") {
    samples <- sampler(n)
    values <- g(samples)
    estimate <- mean(values)
    se <- sd(values) / sqrt(n)

  } else if (method == "antithetic") {
    # Generate n/2 samples and their antithetic pairs
    u <- runif(n/2)
    samples1 <- qnorm(u)
    samples2 <- qnorm(1 - u)  # Antithetic

    values1 <- g(samples1)
    values2 <- g(samples2)
    paired_means <- (values1 + values2) / 2

    estimate <- mean(paired_means)
    se <- sd(paired_means) / sqrt(n/2)

  } else if (method == "control") {
    samples <- sampler(n)
    values <- g(samples)

    # Use sample mean as control (known E[X] = 0 for standard normal)
    control <- samples
    c_star <- -cov(values, control) / var(control)

    adjusted <- values + c_star * (control - 0)
    estimate <- mean(adjusted)
    se <- sd(adjusted) / sqrt(n)
  }

  list(
    estimate = estimate,
    se = se,
    ci = estimate + c(-1.96, 1.96) * se,
    n = n,
    method = method
  )
}

---

Importance Sampling

Theory

To estimate $E_P[g(X)]$ when sampling from $P$ is difficult, sample from proposal $Q$:

$$E_P[g(X)] = E_Q\left[g(X) \frac{p(X)}{q(X)}\right] = E_Q[g(X) w(X)]$$

where $w(X) = p(X)/q(X)$ are importance weights.

Self-Normalized Importance Sampling

When normalizing constants are unknown:

$$\hat{\mu} = \frac{\sum_{i=1}^N w_i g(X_i)}{\sum_{i=1}^N w_i}$$

Effective Sample Size

$$\text{ESS} = \frac{(\sum_i w_i)^2}{\sum_i w_i^2}$$

Rule of thumb: ESS > N/2 indicates reasonable proposal.

R Implementation

#' Importance Sampling Estimator
#'
#' @param g Function to evaluate
#' @param log_target Log of target density (unnormalized OK)
#' @param log_proposal Log of proposal density
#' @param proposal_sampler Function to sample from proposal
#' @param n Number of samples
#' @return Importance sampling estimate
importance_sampling <- function(g, log_target, log_proposal,
                                 proposal_sampler, n = 10000) {

  # Sample from proposal
  samples <- proposal_sampler(n)

  # Compute log importance weights
  log_weights <- log_target(samples) - log_proposal(samples)

  # Stabilize: subtract max for numerical stability
  log_weights <- log_weights - max(log_weights)
  weights <- exp(log_weights)

  # Normalize weights
  normalized_weights <- weights / sum(weights)

  # Compute estimate
  g_values <- g(samples)
  estimate <- sum(normalized_weights * g_values)

  # Effective sample size
  ess <- 1 / sum(normalized_weights^2)

  # Variance estimate (using delta method approximation)
  var_estimate <- sum(normalized_weights^2 * (g_values - estimate)^2)

  list(
    estimate = estimate,
    se = sqrt(var_estimate),
    ess = ess,
    ess_ratio = ess / n,
    max_weight = max(normalized_weights),
    weights = normalized_weights
  )
}

---

MCMC Methods

Metropolis-Hastings Algorithm

**Algorithm**: 1. Initialize $\theta^{(0)}$ 2. For $t = 1, \ldots, T$:

• Propose $\theta^* \sim q(\cdot | \theta^{(t-1)})$
• Compute acceptance probability:

$$\alpha = \min\left(1, \frac{p(\theta^*) q(\theta^{(t-1)} | \theta^*)}{p(\theta^{(t-1)}) q(\theta^* | \theta^{(t-1)})}\right)$$

• Accept with probability $\alpha$

Gibbs Sampling

For multivariate targets, sample each component from its full conditional:

$$\theta_j^{(t)} \sim p(\theta_j | \theta_{-j}^{(t-1)}, \text{data})$$

MCMC Diagnostics

| Diagnostic | Purpose | Target | |------------|---------|--------| | Trace plots | Visual convergence check | Stationary appearance | | $\hat{R}$ (Gelman-Rubin) | Between/within chain variance | < 1.01 | | ESS | Effective independent samples | > 400 per parameter | | Autocorrelation | Mixing assessment | Quick decay |

R Implementation

#' Metropolis-Hastings MCMC
#'
#' @param log_posterior Log posterior function
#' @param init Initial parameter values
#' @param proposal_sd Proposal standard deviation
#' @param n_iter Number of iterations
#' @param n_warmup Warmup iterations to discard
#' @return MCMC samples and diagnostics
metropolis_hastings <- function(log_posterior, init, proposal_sd,
                                 n_iter = 10000, n_warmup = 1000) {

  n_params <- length(init)
  samples <- matrix(NA, nrow = n_iter, ncol = n_params)
  samples[1, ] <- init
  accepted <- 0

  current_lp <- log_posterior(init)

  for (i in 2:n_iter) {
    # Propose
    proposal <- samples[i-1, ] + rnorm(n_params, 0, proposal_sd)