边际似然估计中的数值问题求解

Hey there, let’s work through this tricky numerical stability problem you’re hitting with your Bayesian marginal likelihood estimation!

You’re using an iterative Bayesian framework to estimate marginal likelihood for model selection, where the iteration converges to your target value ( y ). You’ve got vectors ( p_l, q_l ) (length ( L )) and ( p_m, q_m ) (length ( M )), but the core issue is that you’re storing log-likelihood values (with huge negative magnitudes) instead of raw likelihoods. Exponentiating these directly leads to near-zero underflow, and even using a centering constant ( C ) hasn’t fixed things because of the massive scale differences between your vectors.

The key is to rewrite every step of your iteration using log-transformed values, avoiding direct exponentiation entirely. We’ll use the log-sum-exp trick—a standard numerical stability technique—to handle sums of exponential terms without underflow/overflow. Here’s how to adapt your formula:

For each term in your iteration:

1. Convert all raw likelihood references to their stored log-likelihood equivalents.
2. Use log_sum_exp to compute the log of sums of exponential terms (this avoids adding tiny/large numbers directly).
3. Compute the log of each eval_upper/eval_downer term, then use log-sum-exp again to calculate the mean in log space before converting back to a raw value for ( y ).

First, define a robust log_sum_exp function for R:

log_sum_exp <- function(x) {
  max_x <- max(x)
  # Handle cases where all values are -Inf (unlikely for valid likelihoods)
  if (is.infinite(max_x)) return(-Inf)
  max_x + log(sum(exp(x - max_x)))
}

Then update your iteration to use log-space calculations (assuming you’ve loaded your log-likelihood vectors as log_p_l, log_q_l, log_p_m, log_q_m):

L <- 1000 
M <- 5000 
y <- numeric(length = 100)
y[1] <- 0.5  # Initial guess; adjust if needed

for(t in 2:100){
  # Calculate log-space terms for eval_downer and its mean
  log_denom_m <- sapply(1:M, function(m) {
    term1 <- log(L) + log_q_m[m]
    term2 <- log(M) + log_p_m[m] - log(y[t-1])
    log_sum_exp(c(term1, term2))
  })
  log_eval_downer <- log_q_m - log_denom_m
  log_sum_downer <- log_sum_exp(log_eval_downer)
  log_downer <- log_sum_downer - log(M)  # Convert sum log to mean log

  # Calculate log-space terms for eval_upper and its mean
  log_denom_l <- sapply(1:L, function(l) {
    term1 <- log(L) + log_q_l[l]
    term2 <- log(M) + log_p_l[l] - log(y[t-1])
    log_sum_exp(c(term1, term2))
  })
  log_eval_upper <- log_p_l - log_denom_l
  log_sum_upper <- log_sum_exp(log_eval_upper)
  log_upper <- log_sum_upper - log(L)  # Convert sum log to mean log

  # Update y using log-space difference
  y[t] <- exp(log_upper - log_downer)
  cat("Iteration", t, "| Current y:", y[t], "\n")
  
  # Optional: Early stopping if converged
  if (abs(y[t] - y[t-1]) < 1e-8) {
    cat("Converged at iteration", t, "\n")
    break
  }
}

• Early Stopping: Adding a convergence check (like the one above) saves computation time once ( y ) stabilizes.
• Initial Guess: If 0.5 doesn’t work well, you can compute a rough initial ( y ) using log-likelihoods (e.g., a simple harmonic mean estimate in log space).
• Vectorization: For faster computation, replace the sapply calls with vectorized operations (e.g., using purrr::map_dbl or matrix operations) if your dataset is extremely large.

内容的提问来源于stack exchange，提问作者yrx1702