📊 statistics

Here Be Dragons: Bimodal posteriors arise from numerical integration error in longitudinal models

The paper demonstrates that numerical integration errors in longitudinal differential equation models can create artificial bimodality in posterior distributions, leading to the emergence of incorrect parameter estimates that are statistically indistinguishable from the true values.

Original authors: Tess O'Brien, Matthew T. Moores, David Warton, Daniel Falster

Published 2026-02-10

📖 4 min read☕ Coffee break read

CC BY 4.0

Original authors: Tess O'Brien, Matthew T. Moores, David Warton, Daniel Falster

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

The Tale of the Two Maps: Why Computers Sometimes Lie to Us

Imagine you are a treasure hunter trying to find a hidden chest in a dense, foggy forest. You have a map, but the map isn't a drawing; it’s a set of instructions: "Walk 10 steps, turn slightly left, walk 5 steps, repeat."

You follow these instructions carefully. However, there is a catch: you aren't walking in a smooth, continuous glide. You are walking in clunky, robotic jumps. You take one big step, stop, recalculate, and take another.

In this story, the "treasure" is the true biological growth of an animal (like how a fish grows over time), and the "instructions" are mathematical equations. The "robotic jumps" are what scientists call numerical integration—the way a computer approximates a smooth curve by breaking it into tiny, jagged steps.

The Problem: The "Ghost" Treasure

The researchers discovered something spooky. Because the computer is "jumping" through the math instead of sliding smoothly, it can get tricked.

Sometimes, the computer follows the instructions and ends up at the Real Treasure (the true biological growth). But occasionally, because of the way the "jumps" overshoot or undershoot the path, the computer thinks it has found a Ghost Treasure.

This Ghost Treasure is a set of numbers that looks perfectly correct to the computer. If you ask the computer, "Does this look like the right spot?" it says, "Yes! The math checks out!" But in reality, these numbers are total nonsense—they describe a creature that grows at a ridiculous, impossible speed, only to suddenly "teleport" or oscillate wildly to make the math work.

The "Bimodal" Trap: Two Truths, One Lie

In science, we use a method called MCMC to find the answer. Think of MCMC as a team of 10,000 explorers sent into the forest to find the treasure.

Normally, you’d expect all 10,000 explorers to meet up at the same spot (the true answer). But because of this numerical error, the explorers split into two groups:

Group A finds the real treasure.
Group B finds the Ghost Treasure.

When the scientists look at the results, they see two distinct "clusters" of answers. This is what they call bimodality. It’s like a GPS telling you that you are either in London or on the Moon. Both answers are mathematically "valid" according to the computer's broken logic, but only one is real.

Why is this happening? (The Seesaw Metaphor)

The paper uses a specific type of equation (a linear first-order ODE) that acts like a self-correcting seesaw.

If the seesaw tips too far to one side, the math tries to pull it back. But because the computer is taking "robotic jumps," it doesn't just pull the seesaw back to the middle; it slams it to the other side. The computer sees this violent slamming and thinks, "Aha! If I use these crazy, huge numbers, the slamming will perfectly cancel itself out and look like a smooth movement!"

The computer mistakes chaos and error for stability and truth.

How do we fix it?

The researchers suggest three ways to avoid being fooled by the "Dragons" (the errors) in the math:

The "Smooth Slide" (The Analytic Solution): Instead of making the computer take "robotic jumps," give it a formula that allows it to "slide" smoothly along the curve. This bypasses the jumping error entirely.
The "Smaller Steps" Approach: If you must jump, take much, much smaller steps. This makes the jumps so tiny that the "slamming" effect doesn't happen. However, the paper warns that even small steps can still leave you stranded with a Ghost Treasure.
The "Double Check": If you see your explorers splitting into two groups, don't trust the results. Run the math again using a more expensive, more careful method to see if the Ghost Treasure disappears.

The Moral of the Story

Just because a computer says the math "fits" doesn't mean the biology is real. Sometimes, the very tools we use to measure the world create illusions that look exactly like the truth.

Technical Summary: "Here Be Dragons: Bimodal posteriors arise from numerical integration error in longitudinal models"

1. The Problem: Numerical Non-Identifiability

The paper identifies a novel form of non-identifiability in Bayesian inverse methods. Traditionally, non-identifiability arises from insufficient data, large measurement error, or over-parameterization. However, the authors demonstrate that when longitudinal models are governed by differential equations (ODEs) and solved via numerical integration (specifically the Runge-Kutta 4th order (RK4) method), the numerical error itself can create a "phantom" mode in the posterior distribution.

This error manifests as bimodality: the Markov Chain Monte Carlo (MCMC) sampler converges to two distinct sets of parameters. One set represents the true biological/physical parameters, while the second set consists of "erroneous" parameters that are mathematically incorrect but produce numerical projections that are statistically indistinguishable from the true solution due to integration instability.

2. Methodology

The authors use a combination of theoretical proof and empirical simulation to isolate this phenomenon.

Theoretical Framework: They analyze a linear first-order ODE: $\frac{dY}{dt} = \beta_0 - \beta_1 Y(t)$ . This model is chosen because it is a standard component of growth functions (like the von Bertalanffy model) and has a known analytic solution.
Mathematical Proof: Using the RK4 algorithm, the authors derive a polynomial (Proposition 2.1) that describes the relationship between the estimated parameters and the numerical step. They prove that for a given asymptotic size $\alpha$ , there exist multiple values of $\beta_1$ that satisfy the numerical integration step, effectively proving that the numerical method cannot distinguish between the true parameter and the erroneous one.
Empirical Simulation:
- Data Generation: They simulate "true" longitudinal data using the analytic solution with added Gaussian measurement error.
- Estimation: They run 10,000 independent MCMC chains using the hmde framework in Stan to test for bimodality.
- Comparative Testing: They vary the numerical step size ( $h = 0.5, 0.25, 0.125$ ), test different integration methods (RK4 vs. adaptive RK45), and test different prior distributions (Normal vs. Log-normal).
- Validation: They use finite mixture models to statistically confirm the existence of clusters (modes) in the posterior.

3. Key Contributions

Identification of a New Pathology: The paper defines a previously undescribed source of non-uniqueness in parameter estimation: numerical integration error-induced bimodality.
Theoretical Derivation: They provide a formal proof that RK4 integration of linear ODEs inherently allows for multiple parameter solutions that yield the same discrete-time observations.
Simulation Workflow: They propose a qualitative simulation-based workflow that allows researchers to test if their specific model structure is susceptible to these "dragons" (erroneous modes) before committing to large-scale analysis.
Practical Solutions: They demonstrate that using an analytic solution (where possible) or an adaptive step-size method (RK45) can mitigate or eliminate the problem.

4. Results

Persistence of Bimodality: For fixed-step RK4, bimodality is persistent. As the step size decreases, the "bad" mode moves further away from the true value in parameter space, but the probability of the sampler falling into that mode remains significant.
Failure of Priors: Standard techniques, such as adjusting the location or spread of priors, generally fail to eliminate the erroneous mode. The "bad" mode is a structural consequence of the numerical method, not a lack of prior information.
Instability of Erroneous Modes: The erroneous parameter combinations produce highly unstable intermediate sub-steps (oscillating gradients), even though the final projected values at observation times $t_j$ appear correct.
Success of Analytic Solutions: Implementing the model using its analytic solution (avoiding numerical integration entirely) resulted in accurate, unimodal posteriors with negligible bias.

5. Significance

This paper serves as a critical warning for the Bayesian community, particularly in fields like ecology, pharmacology, and physics, where longitudinal models are common. It highlights that numerical methods must be treated as part of the statistical model.

The significance lies in the realization that "well-behaved" MCMC diagnostics (like $\hat{R}$ or trace plots) may fail to detect the issue if independent chains converge to different, yet internally consistent, modes. The authors conclude that researchers must verify the stability of their numerical integrators at the estimated parameter values to ensure that the resulting posterior is a reflection of the underlying biology rather than a mathematical artifact of the solver.