Numerical Stability in Fine-Mapping Calculations

Why Numerical Stability Matters

Fine-mapping calculations can involve likelihood ratios that are extremely small or large. A computer may represent such quantities as floating-point numbers with a limited range and finite precision. Directly exponentiating a large log Bayes factor can overflow; exponentiating a very negative one can underflow to zero. Either problem can corrupt normalized posterior probabilities even when the mathematical calculation is well defined.

The goal here is to preserve the intended calculation while changing how it is implemented. Numerical stabilization does not change the model, the priors, or the evidence. It reduces avoidable computational error in converting those inputs into posterior probabilities and cumulative credible-set mass.

Log Weights and Posterior Normalization

Suppose a finite set of mutually exclusive and exhaustive configurations \(c=1,\ldots,K\) is under consideration. Let \(\pi_c\) be the prior probability of configuration \(c\), let \(BF_c\) be its Bayes factor relative to a common reference model, and assume \(\pi_c\geq 0\), \(BF_c>0\), and \(\sum_{c=1}^K\pi_c=1\). Bayes' rule gives the posterior probability of configuration \(c\) conditional on the data \(D\). The common reference factor cancels in the normalization.

The direct expression can be unstable when \(a_c\) is far outside the range where floating-point exponentiation is reliable. Instead, let \(m=\max_c a_c\), subtract this common value, and exponentiate only the resulting nonpositive differences. Define the log-sum-exp operation as follows.

It follows that the log posterior probability and posterior probability can be computed as

This procedure is stable against overflow in the final normalization because \(a_c-m\leq 0\), so no exponential exceeds 1. Very small terms can still underflow to zero, but only when their contribution is smaller than the computer can represent; they cannot displace the largest term or make the denominator overflow.

Why Subtracting the Maximum Works

The justification is the invariance of a normalized ratio to a common multiplicative factor. Since \(w_c=e^{a_c}\), multiplying every weight by \(e^{-m}\) leaves every ratio unchanged:

Because \(m\) is the largest log weight, \(a_j-m\leq 0\) for all \(j\). Thus each term in the stabilized sum lies between 0 and 1, and at least one term equals 1. The denominator is consequently at least 1, avoiding the particularly dangerous case in which all exponentiated weights vanish. This derivation requires no approximation: only finite-precision arithmetic in the actual implementation introduces numerical error.

Extreme Values: Avoiding Overflow

The same method is especially useful when all log weights are large. A large log weight need not indicate a problem with the posterior: posterior probabilities depend on differences between log weights, not on their common level. Exponentiating the original values may be impossible even though the normalized probabilities are ordinary numbers.

Credible-Set Mass and Accumulation

A credible set is formed from posterior probability assigned to candidate variants or other relevant events. For a single-causal-variant posterior over candidates \(1,\ldots,J\), let \(p_j\) be the posterior probability for candidate \(j\), with \(p_j\geq0\) and \(\sum_j p_j=1\). Sort candidates so that \(p_{(1)}\geq p_{(2)}\geq\cdots\geq p_{(J)}\). For a target level \(q\), such as 0.95, include the fewest highest-probability candidates whose cumulative mass reaches or exceeds \(q\). In settings where candidate probabilities are obtained by summing configuration probabilities, perform that aggregation before ranking candidates.

Use unrounded probabilities for ranking and accumulation. Rounding each term first can change the order of near-ties or move a cumulative sum across the target. For a long list, compensated summation or extended precision can reduce accumulated floating-point error. When probabilities are extremely small, an alternative is to accumulate their log probabilities with the two-term log-sum-exp operation, rather than exponentiating them first.

Practical Checks and Failure Modes

A stable formula cannot rescue incorrect inputs. Before interpreting the output, check that the candidate configurations are the intended mutually exclusive and exhaustive set, the priors are nonnegative, and prior-weighted Bayes factors correspond to the same data and reference model. The log-sum-exp step addresses the numerical normalization; it does not verify those statistical assumptions.

• Overflow in direct exponentiation: intermediate weights become infinite even though their normalized ratios are finite. Subtract the maximum log weight before exponentiating.
• Underflow of small terms: a tiny weight may become zero. If its true contribution could affect a reported result, retain log probabilities and use log-sum-exp, or use higher precision.
• Unstable cancellation in a log Bayes factor: obtaining a log Bayes factor by subtracting two very large, nearly equal log marginal likelihoods can lose significant digits. Use a stable analytic expression or a numerical method designed for that difference, and compare precision when results are consequential.
• Rounding before ranking or summing: rounded values can alter ties and threshold decisions. Keep full precision internally and round only for display.
• Invalid normalization: posterior probabilities should be finite, nonnegative, and sum to 1 within numerical tolerance. A zero or non-finite denominator requires investigation, not an arbitrary rescaling.

For modest log weights, direct exponentiation and normalization may be adequate and easy to inspect. For a wide range of log weights, log-sum-exp is the preferred default. Higher-precision arithmetic can help when many terms or very small probability differences matter, but it is slower and does not correct model misspecification or an inaccurately computed Bayes factor. These approaches differ in numerical behavior, not in the posterior distribution they are meant to represent.