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

TL;DR

This paper reveals that numerical integration error in longitudinal differential-equation models can create a previously undescribed form of non-identifiability, manifesting as bimodal posteriors in Bayesian parameter estimation. The authors prove, for linear first-order ODEs under RK4 integration, that multiple $eta_1$ values can produce indistinguishable trajectories when $eta_0=\alpha\beta_1$, formalized by the root condition $0 = h\hat{\beta}_1 - \frac{h^2}{2}\hat{\beta}_1^2 + \frac{h^3}{6}\hat{\beta}_1^3 - \frac{h^4}{24}\hat{\beta}_1^4 + e^{-\beta_1 h} - 1$, which yields real roots near the true and an erroneous alternative (e.g., approximately 1.0008 and 4.9112 for $\beta_0=10,\beta_1=1$). They corroborate this with empirical simulations across RK4 step sizes, showing two robust posterior modes: one near the true parameter values and another driven by numerical artefacts; switching to an analytic solution eliminates the pathology, while RK45 can mitigate under certain priors. The paper thus provides a practical simulation workflow to detect numerically induced bimodality and warns that numerical integration choices should be treated as part of the statistical model in Bayesian inverse problems. Overall, the work highlights a fundamental interaction between numerical methods and identifiability, with direct implications for longitudinal Bayesian inference and the design of robust estimation workflows.

Abstract

Longitudinal models with dynamics governed by differential equations may require numerical integration alongside parameter estimation. We have identified a situation where the numerical integration introduces error in such a way that it becomes a novel source of non-uniqueness in estimation. We obtain two very different sets of parameters, one of which is a good estimate of the true values and the other a very poor one. The two estimates have forward numerical projections statistically indistinguishable from each other because of numerical error. In such cases, the posterior distribution for parameters is bimodal, with a dominant mode closer to the true parameter value, and a second cluster around the errant value. We demonstrate that multi-modality exists both theoretically and empirically for an affine first order differential equation, that a simulation workflow can test for evidence of the issue more generally, and that Markov Chain Monte Carlo sampling with a suitable solution can avoid bimodality. The issue of multi-modal posteriors arising from numerical error has consequences for Bayesian inverse methods that rely on numerical integration more broadly.

Figures (7)

• Figure 1: Behaviour of numerical solutions at erroneous and true parameter locations showing considerable instability at the erroneous point. The true $Y(t_j)$ values are those of the analytic solution with $\beta_0=10$ and $\beta_1=1$ that the model attempts to estimate, the erroneous values $\beta_0 = 49.3$ and $\beta_1 = 4.91$ can be found both in theoretical results and through empirical testing. (a) and (b) compare the estimated sizes for integer times between 0 and 4 to the analytic solution to show both the good and bad parameters give numerical values very close to the analytic ones. In (c) and (d) we see the sub-steps as well, which are connected up based on the numerical step they come from. The bad parameter estimates demonstrate considerable instability and even a negative value. (e) and (f) give the gradients at the numerical sub-steps and show why the estimated sizes are so unstable for the bad estimate as the gradient oscillates about 0.
• Figure 2: Plot of Equation \ref{['eqn_PropPf3']}, where the $y$-axis replaces 0 on the left hand side, with parameters $h=1/2$, $\beta_1 = 1$, and $\alpha = 10$ for the first five steps. The intersections with 0 are at $\hat{\beta}_1 = 1.0008,\ 4.9112$.
• Figure 3: (a) gives the simulated data with measurement error and precision of 0.1 compared to the analytic solution for $\beta_0 = 10$, $\beta_1 = 1$. Outside of the specified situations with independent error for each chain, estimation is based on the same data. (b) is a scatter plot of all posterior estimates, good and bad, for the default priors across three step sizes using RK4. While there are a total of 30 000 points on this plot, there is so little variation in them that they appear at just four points -- one near the true value, and for each step size, a second smaller cluster of points stacked around an erroneous value. (c) gives the erroneous differential equations and shows that all estimates converge to the same maximum size. (d) shows the analytic solutions using the erroneous parameter combinations to show that the erroneous values give sizes that converge to the asymptotic maximum much faster than the true solution.
• Figure 4: "Well-behaved" chains converging to very different estimates of the parameters $\beta_0 = 10$, $\beta_1 = 1$. (a) produced good estimates while (b) converged to the erroneous parameter combination.
• Figure 5: Estimated $\hat{Y}(t)$ from the linear model across three step sizes compared to the true linear model for simulated data, and numerical error. (a), (c), and (e) show that the estimated sizes over time for the erroneous cluster align closely with the analytic solution for the true $\beta$s across the step sizes. The slight difference in asymptotic size for step size 0.125 is likely due to the prior constraining $\hat{\beta}_0$. (b), (d), and (f) show the numerical error for each erroneous cluster, where the initial step has the largest error and then there is convergence to the same asymptotic size.

Theorems & Definitions (8)

• Lemma 2.1
• proof
• Proposition 2.1
• proof
• Theorem 2.2
• proof
• Corollary 2.1
• proof