R中实现MCMC算法遇接受率为0：泊松-伽马分层模型后验推断排错

Hey there! Let's break down why your MCMC is hitting a 0% acceptance rate—this is a super common snag with hierarchical Gamma-Poisson models, and it almost always traces back to a few key missteps. Here are the most likely culprits to investigate:

If you’re using a random-walk Metropolis-Hastings (MH) algorithm, mismatched proposal scales for your parameters are a top suspect. Since $\lambda$, $\alpha$, and $\beta$ operate on different scales (e.g., $\lambda$ might be a small count rate while $\alpha/\beta$ are hyperparameters with larger values), using the same proposal step size for all three will almost certainly fail.

For example, if you set a large step size for $\alpha$, your candidate values will jump so far from the true posterior that their likelihood is effectively zero, leading to immediate rejection. For positive-only parameters like these, you should also avoid unconstrained normal proposals—they’ll often generate negative candidate values, which have zero posterior density and get rejected automatically.

Fix: Use random walks in the log-space of your parameters (e.g., sample proposals for $\log(\lambda)$, $\log(\alpha)$, $\log(\beta)$) to keep candidates positive, and tune each parameter’s proposal standard deviation independently. Aim for an acceptance rate between 20-50% per parameter.

This is the most frequent mistake in hierarchical models. Your full posterior density should be:
$$p(\lambda, \alpha, \beta | y) \propto p(y|\lambda) \cdot p(\lambda|\alpha,\beta) \cdot p(\alpha) \cdot p(\beta)$$

Common errors here include:

• Forgetting to include the hyperpriors $p(\alpha)$ and $p(\beta)$ in your density calculation
• Misusing R’s dgamma function (remember: R defaults to shape-rate parameterization, not shape-scale—mixing these up will completely break your density)
• Calculating likelihoods as raw products instead of log-sums (this causes numerical underflow, turning valid densities into zero and making acceptance probabilities impossible to compute)

Fix: Write a log-posterior function that sums the log-likelihood, log-prior for $\lambda$, and log-hyperpriors for $\alpha$ and $\beta$. Return -Inf immediately if any parameter is non-positive to avoid invalid calculations.

If you’re sampling in log-space (which you should be for positive parameters), you need to account for the Jacobian determinant when converting back to the original parameter space. For a transformation $x = \log(\theta)$, the Jacobian is $1/\theta$, so your log-posterior must include an additional term of $-\log(\theta)$ (or equivalently, $x$ when working in log-space).

Forgetting this adjustment will make your MH acceptance ratio incorrect, often leading to 0% acceptance because the candidate’s "effective" posterior density is miscalculated.

If you start your MCMC chain far from the bulk of the posterior distribution (e.g., setting an initial $\lambda$ of 100 when the true posterior centers around 5), even small proposal jumps will land in regions with near-zero density. The chain will never move because every candidate is rejected.

Fix: Initialize parameters using reasonable estimates: use the sample mean for $\lambda$, and set $\alpha/\beta$ to values that make the Gamma prior for $\lambda$ match that initial estimate (e.g., $\alpha=2$, $\beta=\lambda_{\text{init}}/2$).

Here’s a simplified, correct log-posterior and MCMC setup to guide you:

# Log-posterior function
log_posterior <- function(params, y) {
  lambda <- params[1]
  alpha <- params[2]
  beta <- params[3]
  
  # Hyperprior parameters (tune these to your knowledge)
  a0 <- 2; b0 <- 1; c0 <- 2; d0 <- 1
  
  # Reject non-positive parameters immediately
  if (lambda <= 0 || alpha <= 0 || beta <= 0) return(-Inf)
  
  # Calculate log-components
  log_lik <- sum(dpois(y, lambda, log = TRUE))
  log_prior_lambda <- dgamma(lambda, shape = alpha, rate = beta, log = TRUE)
  log_prior_alpha <- dgamma(alpha, shape = a0, rate = b0, log = TRUE)
  log_prior_beta <- dgamma(beta, shape = c0, rate = d0, log = TRUE)
  
  return(log_lik + log_prior_lambda + log_prior_alpha + log_prior_beta)
}

# MCMC with log-space sampling
run_mcmc <- function(y, n_iter = 10000, burn_in = 2000) {
  # Initialize with reasonable values
  lambda_init <- mean(y)
  alpha_init <- 2
  beta_init <- lambda_init / alpha_init
  
  # Convert to log-space
  log_lambda <- log(lambda_init)
  log_alpha <- log(alpha_init)
  log_beta <- log(beta_init)
  
  results <- matrix(NA, n_iter, 3)
  colnames(results) <- c("lambda", "alpha", "beta")
  
  # Tunable proposal standard deviations
  prop_sd <- c(0.1, 0.1, 0.1)
  accept_counts <- c(0, 0, 0)
  
  for (i in 1:n_iter) {
    # Update log_lambda with Jacobian adjustment
    log_lambda_cand <- rnorm(1, log_lambda, prop_sd[1])
    lambda_cand <- exp(log_lambda_cand)
    log_p_current <- log_posterior(c(exp(log_lambda), exp(log_alpha), exp(log_beta)), y)
    log_p_cand <- log_posterior(c(lambda_cand, exp(log_alpha), exp(log_beta)), y)
    
    # Jacobian term: log(q(current | candidate)/q(candidate | current)) = log_lambda - log_lambda_cand
    log_accept <- log_p_cand - log_p_current + (log_lambda - log_lambda_cand)
    if (log(runif(1)) < log_accept) {
      log_lambda <- log_lambda_cand
      accept_counts[1] <- accept_counts[1] + 1
    }
    
    # Repeat similar updates for log_alpha and log_beta...
    
    results[i,] <- c(exp(log_lambda), exp(log_alpha), exp(log_beta))
  }
  
  # Print acceptance rates for tuning
  cat("Acceptance Rates:\n")
  cat(paste0("lambda: ", round(accept_counts[1]/n_iter, 3), "\n"))
  cat(paste0("alpha: ", round(accept_counts[2]/n_iter, 3), "\n"))
  cat(paste0("beta: ", round(accept_counts[3]/n_iter, 3), "\n"))
  
  return(results[(burn_in+1):n_iter,])
}

Start by verifying your log-posterior calculation, then tune your proposal scales until you see reasonable acceptance rates. That should get your chain moving!

内容的提问来源于stack exchange，提问作者Amanda R.