REVIEW 3 major objections 6 minor 31 references
Quantifying uncertainty in the numerical integration of evolution equations based on Bayesian isotonic regression
T0 review · 3 major / 6 minor · reviewed 2026-08-12 · deepseek-v4-flash
Pith's one-line read The paper claims that discretization errors in ODE solvers can be quantified as monotone error variances using Bayesian isotonic regression with a horseshoe-type prior, and that approximating the log-chi-square distribution by a…
desk verdict The Gibbs sampler as printed does not target the posterior it claims; the paper's core algorithm is wrong as written, but the modeling idea is solid and worth a serious referee. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
The engine is the variable transformation $\eta_1 = \log\sigma_1^2$ and $\eta_j = \log\sigma_j^2 - \log\sigma_{j-1}^2$, which turns the isotonic constraint $\gamma^2 \le \sigma_1^2 \le \cdots \le \sigma_n^2$ into lower-tail constraints on real-valued $\eta_j$'s. On these increments the paper puts a horseshoe-type scale-mixture prior: $\eta_1 \mid \tau_1 \sim N_{\ge\log\gamma^2}(\log\gamma^2,\tau_1)$ and $\eta_j \mid \tau_j,\lambda \sim N_{\ge 0}(0,\lambda\tau_j)$, with half-Cauchy marginals for $\sqrt{\tau_j}$ and $\sqrt{\lambda}$ giving strong shrinkage with occasional jumps. The second mechanism is the log-square trick: writing $z_i = \log r_i^2$ makes the cumulative $\eta$'s appear as location parameters in a regression of $z_i$ on log-$\chi^2_1$ noise. Replacing that noise density by a 10-component Gaussian mixture with tabulated weights, means, and variances lets the sampler treat the problem as conditionally Gaussian, so all Gibbs full conditionals become truncated normal, gamma, or generalized inverse Gaussian draws. This yields an $O(n^2)$-per-iteration sampling algorithm for the posterior of the discretization error variances.
What would settle it
Compute the actual discretization error sequence for the simple harmonic oscillator integrated with the Störmer–Verlet method over many periods, where the error variance oscillates rather than increasing monotonically. If the proposed method's 90% credible intervals cover the true error in substantially fewer than 90% of time points, the monotone-variance model is contradicted.
Extended reading notes
Core claim
The central claim is that discretization error variances $\sigma_i^2$ of a numerical ODE integrator can be estimated from observable residuals $r_i = v_i - O(x_i)$ by a Bayesian isotonic regression model, provided the variances are assumed nondecreasing and the observation noise variance $\gamma^2$ is known. Writing $\eta_1 = \log\sigma_1^2$ and $\eta_j = \log\sigma_j^2 - \log\sigma_{j-1}^2$ for $j\ge 2$ converts monotonicity into the simple constraints $\eta_1 \ge \log\gamma^2$ and $\eta_j \ge 0$. Placing a half-shrinkage (horseshoe) prior on the $\eta_j$'s induces strong shrinkage of adjacent variances toward equality while still allowing occasional sharp jumps. The log-squared residuals satisfy $z_i = \log r_i^2 = \eta_1 + \cdots + \eta_i + \log\varepsilon_i^2$, where $\log\varepsilon_i^2$ has the log-$\chi^2_1$ distribution; approximating that density by a 10-component Gaussian mixture with fixed weights, means, and variances makes every full conditional in the Gibbs sampler a standard truncated normal, gamma, or generalized inverse Gaussian distribution. The paper demonstrates on the FitzHugh–Nagumo model and a Kepler-orbit velocity observation that the resulting posterior means and 95% credible intervals for $\sigma_i$ and for $\sqrt{\sigma_i^2-\gamma^2}$ track the actual residual and error scales.
Load-bearing premise
The load-bearing premise is that discretization errors are independent, zero-mean Gaussian random variables whose variances never decrease; the paper itself flags this in Remark 2, and if the true error is correlated, biased, or has a decreasing phase, the posterior credible intervals will be miscalibrated.
Editorial extensions
If this is right
- Posterior means and 95% credible intervals for each $\sigma_i$ (and for the pure discretization part $\sqrt{\sigma_i^2-\gamma^2}$) provide an uncertainty band around the numerical solution at every observation time.
- Because all full conditionals are standard distributions, the Gibbs sampler is straightforward to implement, making error quantification practical when the ODE state dimension is large relative to the number of observations.
- The construction extends Bayesian isotonic regression from estimating normal means to estimating normal variances, so it can be reused for other problems with a monotone variance structure.
- In the FitzHugh–Nagumo and Kepler experiments, a large proportion of observed residuals and true errors fall inside the 90% credible intervals, indicating the estimated variances capture the error scale.
Reading between the lines
- A natural next step is to couple this error quantification with Bayesian parameter estimation for the ODE itself, so that discretization uncertainty is propagated into posterior distributions of unknown parameters rather than being ignored.
- The diagonal-covariance restriction could be lifted using a Wishart-type prior on full matrices, an extension the paper explicitly leaves for future work.
- The fixed 10-component mixture approximation could be audited by comparing posterior intervals under more components or alternative approximations, and the posterior output could serve as a diagnostic for whether the monotone-variance assumption actually holds.
- Because the method yields a full posterior over error variances, it could be used to choose step sizes adaptively by shrinking steps where the estimated variance is large, complementing traditional local error control.
Signed reviews
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The manuscript proposes a Bayesian framework for quantifying discretization errors in numerical solutions of ordinary differential equations. The residuals between noisy observations and numerical approximations are modeled as independent zero-mean Gaussians with monotonically nondecreasing variances, and a horseshoe-type shrinkage prior is placed on the log-increments of these variances. To make posterior sampling tractable, the log-chi-square distribution arising from a log-squared transformation is approximated by a 10-component Gaussian mixture, and a Gibbs sampler is constructed. The method is demonstrated on the FitzHugh-Nagumo and Kepler models, where posterior credible intervals are compared with actual residual and error values.
Significance. The paper addresses an important problem in probabilistic numerics: giving uncertainty quantification for deterministic ODE discretization errors. Extending Bayesian isotonic regression from normal means to variances is a useful conceptual step, and the Gaussian-mixture approximation of the log-chi-square distribution is a sensible device for enabling conjugate-style updates. If the sampling algorithm is corrected, the framework could be practically valuable for inverse problems in ODE models. The authors are also explicit about the strong modeling assumptions (Gaussianity, independence, monotone variances), which is commendable. However, the central algorithmic claim is undermined by an incorrect full-conditional update for the mixture indicators, so the numerical results as presented cannot be used to validate the stated posterior.
major comments (3)
- [Section 2.3.2 and Algorithm 1] The update for the auxiliary mixture indicators s_i is not a valid Gibbs step. The paper states that s_i can be sampled from the prior probability mass (5), and Algorithm 1 step 4 draws s(t) ~ P(s). In the augmented model, the full conditional is p(s_i = k | z_i, eta_i) proportional to w_k N(z_i - eta_i; m_k, v_k^2), which depends on the data through z_i - eta_i. Because the subsequent update of eta_j in step 5-6 uses the current values of tau and lambda, the overall transition kernel depends on the current state in a way that does not have the intended posterior as its stationary distribution. Thus the sampler as printed does not target pi(Sigma, Upsilon | r), and the coverage results in Section 3 cannot be attributed to the stated posterior.
- [Section 2.3.2, formula for w^2_j] The variance in the conditional posterior of eta_j is printed as w^2_j = (tilde{w}^2_j lambda tau_1) / (tilde{w}^2_j + lambda tau_j), and the same expression is repeated in Algorithm 1 line 10. Since the prior for eta_j is N_{>=0}(0, lambda tau_j), the prior precision is 1/(lambda tau_j), so the correct posterior precision is 1/tilde{w}^2_j + 1/(lambda tau_j), giving w^2_j = (tilde{w}^2_j lambda tau_j) / (tilde{w}^2_j + lambda tau_j). The numerator should contain tau_j, not tau_1. This is a concrete formula error that would bias the sampling of eta_j if the algorithm is implemented as written.
- [Section 3, coverage assessment] The coverage checks in Figures 5-8 are in-sample: the same residuals are used to estimate the variance parameters and then to check whether the 90% credible intervals contain the residuals or errors. The text reports only that 'a large proportion of data points fall within the 90% credible intervals' for single trajectories, without a quantitative coverage rate. Because the model is fit to the same data being evaluated, the intervals are expected to show high coverage by construction. To support the central claim that the method accurately quantifies discretization error, the authors should report repeated-simulation coverage (e.g., the fraction of held-out errors covered across many trajectories) or otherwise assess out-of-sample behavior.
minor comments (6)
- [Section 2.3.2] The sentence 'Since the latent parameters s are independent of other parameters' is only true in the prior; in the posterior, s is correlated with eta and the data. This wording facilitates the incorrect update described in the first major comment.
- [Algorithm 1, line 10] The line reads 'Sample tau^{(t)}_1 ~ P(tau_j | eta^{(t)}_j, nu^{(t)}_j)'; the subscript on the sampled variable should be j, not 1.
- [Figure 1 caption] The caption contains a typo: 'the Gaussian mixture model g(epsilon) give in (4)' should read 'given in (4)'.
- [Remark 4] The phrase 'In such a case case' contains a duplicated word and should be 'In such a case'.
- [Section 2.3.1] The GIG parameter in the update for lambda is typeset as '-n + 2/2', which is ambiguous; it should be '-(n-2)/2'.
- [Section 3] No convergence diagnostics (trace plots, R-hat) are reported for the Gibbs sampler, and no code is released for the proposed method. Given the algorithmic error in the printed sampler, reproducible code and convergence checks would be essential for validation.
Circularity Check
No significant circularity: the Bayesian model is built from explicit priors and standard mixture approximations, and the in-sample coverage checks are model calibration, not disguised predictions.
full rationale
The central derivation chain is a Bayesian hierarchical model: residuals r_i = v_i - O(x_i) are modeled as N(0, sigma_i^2) with sigma_i^2 = gamma^2 + V_i (Section 2.1); the monotone variances are parameterized by log-increments eta_j with a horseshoe-type shrinkage prior adapted from Okano et al. [16]; and the log-chi-square density is replaced by the published 10-component Gaussian mixture of Omori et al. [29] (Section 2.3.2). These are explicit statistical inputs (prior and approximation), not outputs of the inference, so the posterior is defined by Bayes theorem rather than being equivalent to its inputs by construction. The numerical checks in Section 3 plot posterior credible intervals for the absolute residual |r_i| and for the discretization error sqrt(sigma_i^2 - gamma^2), then observe that a large proportion of the data points fall inside. These are in-sample posterior predictive checks, not out-of-sample predictions, so the intervals are not 'predictions' in the sense of being statistically forced by construction; moreover, the discretization error itself (x(t_i)-x_i) is not directly fed into the likelihood, giving the error check partial diagnostic content. Self-citations [12-14] and [16] (with co-author Irie) are used for background, motivation, and prior choice, but they are not load-bearing justifications for the central claim; Remark 2 explicitly acknowledges the Gaussian zero-mean assumption as a simplification. Two non-circular correctness issues should be noted: Section 2.3.2 and Algorithm 1 sample s_i from the prior P(s) rather than the conditional posterior p(s_i|z,eta), and the printed variance formula for eta_j has tau_1 in the numerator where tau_j is required. These may affect the numerical results but are not circularity. Overall, no step reduces by construction to its input, so the circularity score is low.
Assumptions & free parameters
free parameters (2)
- sigma_i^2 (i=1,...,n) =
posterior means displayed in Figures 3-8
- Prior hyperparameters (shapes and rates of Gamma priors on nu_j and xi) =
fixed at 1/2 or 1
assumptions (5)
- domain assumption Discretization errors O(x(t_i))-O(x_i) are independent zero-mean Gaussian with variance V_i.
- domain assumption The variances satisfy 0 < gamma^2 <= sigma_1^2 <= sigma_2^2 <= ... <= sigma_n^2.
- domain assumption Observation noise covariance Gamma is known and diagonal, and the observation operator is handled directly.
- domain assumption The 10-component Gaussian mixture g(epsilon_tilde) accurately approximates the log-chi-square_1 density f(epsilon_tilde).
- standard math The Gibbs sampler converges to the target posterior, and GIG variate generation is reliable in the edge cases discussed in Remark 4.
Cite this review
Pith. "Pith review of Quantifying uncertainty in the numerical integration of evolution equations based on Bayesian isotonic regression." pith.science (2026). https://pith.science/paper/PGN33AOH
@misc{pith2026241108338,
author = {Pith},
title = {Pith review of: Quantifying uncertainty in the numerical integration of evolution equations based on Bayesian isotonic regression},
year = {2026},
howpublished = {\url{https://pith.science/paper/PGN33AOH}},
note = {Machine review of arXiv:2411.08338}
}
abstract
This paper presents a new Bayesian framework for quantifying discretization errors in numerical solutions of ordinary differential equations. By modelling the errors as random variables, we impose a monotonicity constraint on the variances, referred to as discretization error variances. The key to our approach is the use of a shrinkage prior for the variances coupled with variable transformations. This methodology extends existing Bayesian isotonic regression techniques to tackle the challenge of estimating the variances of a normal distribution. An additional key feature is the use of a Gaussian mixture model for the $\log$-$\chi^2_1$ distribution, enabling the development of an efficient Gibbs sampling algorithm for the corresponding posterior.
Reference graph
Works this paper leans on
-
[1]
II: Stiff and Differential-Algebraic Problems
Hairer, E., Wanner, G.: Solving Ordinary Differential Equations. II: Stiff and Differential-Algebraic Problems. Springer Series in Computational Mathematics, vol. 14. Springer, Berlin (1991). https://doi.org/10.1007/978-3-662-09947-6
-
[2]
I: 16 Nonstiff Problems, 2nd edn
Hairer, E., Nørsett, S.P., Wanner, G.: Solving Ordinary Differential Equations. I: 16 Nonstiff Problems, 2nd edn. Springer Series in Computational Mathematics, vol
-
[3]
Conrad, P.R., Girolami, M., S¨ arkk¨ a, S., Stuart, A., Zygalakis, K.: Statistical analysis of differential equations: introducing probability measures on numeri- cal solutions. Stat. Comput. 27(4), 1065–1082 (2017) https://doi.org/10.1007/ s11222-016-9671-0
work page 2017
-
[4]
Kersting, H., Sullivan, T.J., Hennig, P.: Convergence rates of Gaussian ODE filters. Stat. Comput. 30(6), 1791–1816 (2020) https://doi.org/10.1007/ s11222-020-09972-4
work page 2020
-
[5]
Schober, M., S¨ arkk¨ a, S., Hennig, P.: A probabilistic model for the numerical solution of initial value problems. Stat. Comput. 29(1), 99–122 (2019) https: //doi.org/10.1007/s11222-017-9798-7
-
[6]
Tronarp, F., Kersting, H., S¨ arkk¨ a, S., Hennig, P.: Probabilistic solutions to ordi- nary differential equations as nonlinear Bayesian filtering: a new perspective. Stat. Comput. 29(6), 1297–1315 (2019) https://doi.org/10.1007/s11222-019-09900-1
-
[7]
Tronarp, F., S¨ arkk¨ a, S., Hennig, P.: Bayesian ODE solvers: the maximum a posteriori estimate. Stat. Comput. 31(3), 23–18 (2021) https://doi.org/10.1007/ s11222-021-09993-7
work page 2021
-
[8]
Springer, Berlin (1993)
work page 1993
Show all 31 references
-
[9]
Abdulle, A., Garegnani, G.: Random time step probabilistic methods for uncer- tainty quantification in chaotic and geometric numerical integration. Stat. Comput. 30(4), 907–932 (2020) https://doi.org/10.1007/s11222-020-09926-w
2020 doi
-
[10]
Lie, H.C., Stuart, A.M., Sullivan, T.J.: Strong convergence rates of probabilistic integrators for ordinary differential equations. Stat. Comput. 29(6), 1265–1283 (2019) https://doi.org/10.1007/s11222-019-09898-6
2019 doi
-
[11]
Calcolo 59(1), 13–33 (2022) https://doi.org/10.1007/s10092-022-00457-6
Lie, H.C., Stahn, M., Sullivan, T.J.: Randomised one-step time integration meth- ods for deterministic operator differential equations. Calcolo 59(1), 13–33 (2022) https://doi.org/10.1007/s10092-022-00457-6
2022 doi
-
[12]
Cambridge University Press, Cambridge (2022)
Hennig, P., Osborne, M.A., Kersting, H.P.: Probabilistic Numerics—Computation as Machine Learning. Cambridge University Press, Cambridge (2022). https:// doi.org/10.1017/9781316681411
2022 doi
-
[13]
Marumo, N., Matsuda, T., Miyatake, Y.: Modelling the discretization error of initial value problems using the Wishart distribution. Appl. Math. Lett. 147, 108833–9 (2024) https://doi.org/10.1016/j.aml.2023.108833
2024
-
[14]
SIAM/ASA J
Matsuda, T., Miyatake, Y.: Estimation of ordinary differential equation models with discretization error quantification. SIAM/ASA J. Uncertain. Quantif. 9(1), 302–331 (2021) https://doi.org/10.1137/19M1278405 17
2021 doi
-
[15]
Matsuda, T., Miyatake, Y.: Piecewise monotone estimation in one-parameter exponential family (2023)
2023
-
[16]
Oates, C.J., Cockayne, J., Aykroyd, R.G., Girolami, M.: Bayesian probabilistic numerical methods in time-dependent state estimation for industrial hydrocyclone equipment. J. Amer. Statist. Assoc. 114(528), 1518–1531 (2019) https://doi.org/ 10.1080/01621459.2019.1574583
2019
-
[17]
Okano, R., Hamura, Y., Irie, K., Sugasawa, S.: Locally adaptive Bayesian isotonic regression using half shrinkage priors. Scand. J. Statist. 51(1), 109–141 (2024) https://doi.org/10.1111/sjos.12676
2024 doi
-
[18]
Cai, B., Dunson, D.B.: Bayesian multivariate isotonic regression splines: Appli- cations to carcinogenicity studies. J. Amer. Statist. Assoc. 102(480), 1158–1171 (2007) https://doi.org/10.1198/016214506000000942
2007 doi
-
[19]
Curtis, S.M., Ghosh, S.K.: A variable selection approach to monotonic regression with Bernstein polynomials. J. Appl. Stat. 38(5), 961–976 (2011) https://doi. org/10.1080/02664761003692423
2011 doi
-
[20]
Bio- metrics 60(2), 398–406 (2004) https://doi.org/10.1111/j.0006-341X.2004.00184
Neelon, B., Dunson, D.B.: Bayesian isotonic regression and trend analysis. Bio- metrics 60(2), 398–406 (2004) https://doi.org/10.1111/j.0006-341X.2004.00184. x
2004 arXiv
-
[21]
Shively, T.S., Sager, T.W., Walker, S.G.: A Bayesian approach to non-parametric monotone function estimation. J. R. Stat. Soc. Ser. B Stat. Methodol. 71(1), 159–175 (2009) https://doi.org/10.1111/j.1467-9868.2008.00677.x
2009
-
[22]
In: Artificial Intelligence and Statistics, pp
Carvalho, C.M., Polson, N.G., Scott, J.G.: Handling sparsity via the horseshoe. In: Artificial Intelligence and Statistics, pp. 73–80 (2009). PMLR
2009
-
[23]
Biometrika 97(2), 465–480 (2010) https://doi.org/10.1093/biomet/asq017
Carvalho, C.M., Polson, N.G., Scott, J.G.: The horseshoe estimator for sparse sig- nals. Biometrika 97(2), 465–480 (2010) https://doi.org/10.1093/biomet/asq017
2010 doi
-
[24]
Dagpunar, J.S.: An easily implemented generalised inverse Gaussian generator. Comm. Statist. Simulation Comput. 18(2), 703–710 (1989) https://doi.org/10. 1080/03610918908812785
1989
-
[25]
Diploma Thesis, Technical University Graz, Austria (1989)
Lehner, K.: Erzeugung von Zufallszahlen aus zwei exotischen Verteilungen. Diploma Thesis, Technical University Graz, Austria (1989)
1989
-
[26]
H¨ ormann, W., Leydold, J.: Generating generalized inverse Gaussian ran- dom variates. Stat. Comput. 24(4), 547–557 (2014) https://doi.org/10.1007/ s11222-013-9387-3
2014
-
[27]
PhD thesis, Massachusetts Institute of Technology (1988) 18
Nelson, D.B.: The time series behavior of stock market volatility and returns. PhD thesis, Massachusetts Institute of Technology (1988) 18
1988
-
[28]
Harvey, A., Ruiz, E., Shephard, N.: Multivariate stochastic variance models. Rev. Econ. Stud. 61(2), 247–264 (1994) https://doi.org/10.2307/2297980
1994 doi
-
[29]
Kim, S., Shephard, N., Chib, S.: Stochastic volatility: likelihood inference and comparison with ARCH Models. Rev. Econ. Stud. 65(3), 361–393 (1998) https: //doi.org/10.1111/1467-937X.00050
1998 doi
-
[30]
Omori, Y., Chib, S., Shephard, N., Nakajima, J.: Stochastic volatility with lever- age: Fast and efficient likelihood inference. J. Econom. 140(2), 425–449 (2007) https://doi.org/10.1016/j.jeconom.2006.07.008
2007 doi
-
[31]
arXiv preprint arXiv:2404.13986 (2024) 19
Hiraki, D., Chib, S., Omori, Y.: Stochastic volatility in mean: Efficient analysis by a generalized mixture sampler. arXiv preprint arXiv:2404.13986 (2024) 19
2024 arXiv
Reviewed August 12, 2026 · model on record in the stance chip above.
Discussion (0). Continue with ORCID to comment.