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

  
  normalized_weights  weights  weights

  
  g_values  gsamples
  estimate  normalized_weights  g_values

  
  ess  1  normalized_weights

  
  var_estimate  normalized_weights  g_values  estimate

  
    estimate  estimate
    se  var_estimate
    ess  ess
    ess_ratio  ess  n
    max_weight  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)

   i  n_iter 
    
    proposal  samplesi   rnormn_params  proposal_sd

    
    proposed_lp  log_posteriorproposal
    log_alpha  proposed_lp  current_lp

    
     runif  log_alpha 
      samplesi   proposal
      current_lp  proposed_lp
       i  n_warmup accepted  accepted  
      
      samplesi   samplesi 
    
  

  
  samples  samplesn_warmup  n_iter 

  
  ess  applysamples  x 
    acf_vals  acfx plot   lag.max  acf
    n  x
    tau  1    acf_vals
    n  tau
  

  
    samples  samples
    acceptance_rate  accepted  n_iter  n_warmup
    ess  ess
    summary  data.frame
      mean  colMeanssamples
      sd  applysamples  sd
      q025  applysamples  quantile 
      q975  applysamples  quantile 
      ess  ess
    
  

---

Hamiltonian Monte Carlo

Theory

HMC uses Hamiltonian dynamics to propose distant moves with high acceptance:

$$H(\theta, p) = -\log p(\theta | y) + \frac{1}{2}p^T M^{-1} p$$

Leapfrog integrator:

1. $p \leftarrow p - \frac{\epsilon}{2} \nabla_\theta U(\theta)$
2. $\theta \leftarrow \theta + \epsilon M^{-1} p$
3. $p \leftarrow p - \frac{\epsilon}{2} \nabla_\theta U(\theta)$

Key Parameters

Parameter	Description	Tuning
$\epsilon$	Step size	Target ~65% acceptance
$L$	Number of leapfrog steps	Balance ESS vs. computation
$M$	Mass matrix	Approximate posterior covariance

Stan Integration

#' Fit Bayesian Mediation Model with Stan
#'
#' @param data List with y, m, x, covariates
#' @return Stan fit object
fit_mediation_stan <- function(data) {
  stan_code <- "
  data {
    int<lower=0> N;
    vector[N] y;
    vector[N] m;
    vector[N] x;
  }
  parameters {
    real a;           // X -> M path
    real b;           // M -> Y path
    real c_prime;     // X -> Y direct
    real alpha_m;     // M intercept
    real alpha_y;     // Y intercept
    real<lower=0> sigma_m;
    real<lower=0> sigma_y;
  }
  model {
    // Priors
    a ~ normal(0, 1);
    b ~ normal(0, 1);
    c_prime ~ normal(0, 1);

    // Likelihoods
    m ~ normal(alpha_m + a * x, sigma_m);
    y ~ normal(alpha_y + b * m + c_prime * x, sigma_y);
  }
  generated quantities {
    real indirect = a * b;
    real total = a * b + c_prime;
  }
  "

  stan_fit <- rstan::stan(
    model_code = stan_code,
    data = data,
    chains = 4,
    iter = 2000,
    warmup = 1000,
    cores = parallel::detectCores()
  )

  stan_fit
}

---

Parallel Computing

Embarrassingly Parallel Tasks

Bootstrap, cross-validation, and simulation studies parallelize easily:

#' Parallel Bootstrap for Mediation
#'
#' @param data Dataset
#' @param statistic Function computing indirect effect
#' @param R Number of bootstrap replicates
#' @return Bootstrap distribution
parallel_bootstrap <- function(data, statistic, R = 2000) {
  library(parallel)

  n <- nrow(data)
  n_cores <- detectCores() - 1

  # Create cluster
  cl <- makeCluster(n_cores)
  clusterExport(cl, c("data", "statistic", "n"), envir = environment())

  # Parallel bootstrap
  boot_results <- parSapply(cl, R i 
    boot_idx  samplen replace  
    boot_data  databoot_idx 
    statisticboot_data
  

  stopClustercl

  
    estimate  statisticdata
    boot_se  sdboot_results
    boot_ci  quantileboot_results  
    boot_dist  boot_results
  

Future Package for Flexible Parallelism

#' Parallel Simulation Study with future
#'
#' @param dgp Data generating process function
#' @param estimator Estimation function
#' @param n_sims Number of simulations
#' @param params Parameter grid
#' @return Simulation results
parallel_simulation <- function(dgp, estimator, n_sims = 1000, params) {
  library(future)
  library(future.apply)

  # Set up parallel backend
  plan(multisession, workers = parallel::detectCores() - 1)

  results <- future_lapply(1:n_sims, function(sim) {
    data <- dgp(params)
    est <- estimator(data
    sim  sim est
   future.seed  

  plansequential  

  do.callrbind results

---

GPU Acceleration

When to Use GPU

Task	CPU	GPU	Recommendation
Small matrices (<1000)	Fast	Overhead	CPU
Large matrices (>5000)	Slow	Fast	GPU
Element-wise operations	Moderate	Very fast	GPU if large
MCMC (sequential)	Fast	Overhead	CPU
Parallel MCMC chains	Moderate	Fast	GPU

R GPU Computing with torch

#' GPU-Accelerated Matrix Operations
#'
#' @param X Design matrix
#' @param y Response vector
#' @return OLS coefficients computed on GPU
gpu_ols <- function(X, y) {
  library(torch)

  # Move to GPU
  X_gpu <- torch_tensor(X, device = "cuda")
  y_gpu <- torch_tensor(y, device = "cuda")

  # Solve normal equations on GPU
  XtX <- torch_mm(torch_t(X_gpu), X_gpu)
  Xty <- torch_mm(torch_t(X_gpu), y_gpu)

  # Solve system
  beta <- torch_linalg_solve(XtX, Xty)

  # Move back to CPU
  as.numeric(beta$cpu(

---

Approximate Bayesian Computation

ABC Rejection Algorithm

When likelihood is intractable but simulation is possible:

1. Sample $\theta^* \sim \pi(\theta)$ (prior)
2. Simulate $y^* \sim p(y | \theta^*)$
3. Accept if $d(S(y^*), S(y_{obs})) < \epsilon$

ABC for Mediation

#' ABC for Complex Mediation Model
#'
#' @param y_obs Observed outcome
#' @param m_obs Observed mediator
#' @param x_obs Observed treatment
#' @param prior_sampler Function to sample from prior
#' @param simulator Function to simulate data given parameters
#' @param summary_stats Function to compute summary statistics
#' @param epsilon Acceptance threshold
#' @param n_samples Target number of accepted samples
#' @return ABC posterior samples
abc_mediation <- function(y_obs, m_obs, x_obs,
                           prior_sampler, simulator, summary_stats,
                           epsilon = 0.1, n_samples = 1000) {

  obs_stats <- summary_stats(y_obs, m_obs, x_obs)
  accepted <- list(
  n_tried  0

   accepted  n_samples 
    
    theta  prior_sampler

    
    sim_data  simulatortheta y_obs
    sim_stats  summary_statssim_datay sim_datam x_obs

    
    distance  sim_stats  obs_stats
     distance  epsilon 
      acceptedaccepted    theta
    

    n_tried  n_tried  
  

  
    samples  do.callrbind accepted
    acceptance_rate  n_samples  n_tried
    epsilon  epsilon
  

---

Convergence Diagnostics

Diagnostic Checklist

• Multiple chains started from dispersed initial values
• $\hat{R} < 1.01$ for all parameters
• ESS > 400 for all parameters
• Trace plots show stationarity
• No divergent transitions (HMC)
• Energy diagnostics pass (HMC)

R Diagnostic Functions

#' Comprehensive MCMC Diagnostics
#'
#' @param samples Matrix of MCMC samples (iterations x parameters)
#' @param n_chains Number of chains (samples split equally)
#' @return Diagnostic summary
mcmc_diagnostics <- function(samples, n_chains = 4) {

  n_iter <- nrow(samples) / n_chains
  n_params <- ncol(samples)

  # Split into chains
  chains <- lapply(1:n_chains, function(i) {
    idx <- ((i-1) * n_iter + 1):(i * n_iter)
    samples[idx, , drop = FALSE]
  

  
  compute_rhat  param_idx 
    chain_means  sapplychains  mean param_idx
    chain_vars  sapplychains  var param_idx

    W  meanchain_vars
    B  varchain_means  n_iter

    var_hat  n_iter    n_iter  W  B  n_iter
    var_hat  W
  

  rhat  sapplyn_params compute_rhat

  
  compute_ess  x 
    n  x
    acf_vals  acfx plot   lag.max  n acf
    tau  1    acf_vals
     n  tau
  

  ess  applysamples  compute_ess

  
    rhat  rhat
    ess  ess
    all_converged  rhat    ess  
    summary  data.frame
      parameter  n_params
      rhat  rhat
      ess  ess
      converged  rhat    ess  
    
  

---

References

Monte Carlo Methods

• Robert, C. P., & Casella, G. (2004). Monte Carlo Statistical Methods
• Owen, A. B. (2013). Monte Carlo theory, methods and examples

MCMC

• Brooks, S., et al. (2011). Handbook of Markov Chain Monte Carlo
• Gelman, A., et al. (2013). Bayesian Data Analysis (3rd ed.)

Computational Statistics

• Gentle, J. E. (2009). Computational Statistics
• Givens, G. H., & Hoeting, J. A. (2012). Computational Statistics

Software

• Stan Development Team. Stan User's Guide
• R Core Team. parallel package documentation

---

Version: 1.0.0 Created: 2025-12-09 Domain: Computational methods for statistical inference Prerequisites: Probability theory, Bayesian statistics, R programming