{"id":"fba87e2f-6881-4c54-a5ac-7d07445274bf","arxiv_id":"2411.08338","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":2,"one_line_summary":"A Bayesian isotonic regression model with a shrinkage prior estimates nondecreasing discretization error variances for numerical ODE solutions using a Gibbs sampler.","lead":"This paper builds a Bayesian method to estimate how wrong a numerical ODE solution is at each time step. It assumes the error grows over time and uses a statistical model to put uncertainty bounds on that error.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The printed Gibbs update for the mixture indicators s_i samples from the prior (5), not the conditional posterior, so as written the sampler does not target π(Σ,Υ|r).","rationale":"The reader's weakest assumption concerns the Gaussian/monotone modeling of discretization errors; that is a real limitation and is explicitly acknowledged in Remark 2. But the more immediate, load-bearing issue for the central algorithmic claim is that the printed Gibbs sampler may not sample from the posterior of even that model. The s update as described is the prior rather than the full conditional; neglecting the likelihood term in the auxiliary mixture step breaks the data augmentation. In addition, the τ_1/τ_j typo in the η_j variance formula would distort the scale of the proposed draws. These are internal inconsistencies, not matters of modeling consensus, and they are checkable. Because the paper ships no code and the numerical section gives no convergence diagnostics or comparison with a reference sampler, the current text cannot establish the central claim. The verdict should remain conditional: accept only after the s update and variance formula are corrected and validated against an independent sampler, or code is released demonstrating the correct implementation.","tokens_in":12377,"tokens_out":17556,"duration_ms":179516,"concrete_test":"Implement Algorithm 1 exactly as printed, including drawing s_i from (5) without conditioning on z and using the printed τ_1 term in the η_j variance, on a small synthetic dataset with known monotonically increasing σ_i^2. Compare the posterior distribution of (σ_1^2,...,σ_n^2) with draws from an independent reference sampler targeting the exact posterior of the same model (e.g., HMC in Stan or random-walk Metropolis on the log-χ^2 likelihood). If the two posterior distributions differ beyond Monte Carlo error, the printed algorithm is not a valid Gibbs sampler. Then change only the s update to p(s_i=k|·) ∝ w_k N(z_i-η_i; m_k, v_k^2) and repeat; matching the reference would confirm the s update is the cause.","verdict_should_be":"UNCHANGED","load_bearing_attack":"Section 2.3.2 states the full conditional for the auxiliary mixture indicators s_i as: 'Since the latent parameters s are independent of other parameters. s_i can be sampled using the uniform distribution with the probability mass (5).' In the augmented model (4)-(5), however, the conditional posterior is p(s_i = k | z_i, η_i) ∝ w_k N(z_i - η_i; m_k, v_k^2); it depends on the observed transformed residual through z_i - η_i. Algorithm 1 step 4 also samples s from P(s) with no dependence on z. A Gibbs sampler that draws s from the marginal prior does not have π(Σ, Υ | r) as its stationary distribution, because the mixture indicators are not conditioned on the data. The same section also contains a concrete formula error: in the update for η_j the printed variance is w^2_j = \\tilde{w}^2_j λ τ_1 / (\\tilde{w}^2_j + λ τ_j), where the numerator should contain τ_j, not τ_1 (the full conditional precision from the prior is 1/(λτ_j)); Algorithm 1 line 10 repeats this τ_1/τ_j slip. Since no code is released, the numerical results cannot show that the sampler as printed converges to the intended posterior.","agreement_with_reader":"disagree"},"referee_report":{"model":"deepseek-v4-flash","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.","tokens_in":12621,"tokens_out":19060,"duration_ms":161471,"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":[{"comment":"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":"Section 2.3.2 and Algorithm 1"},{"comment":"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":"Section 2.3.2, formula for w^2_j"},{"comment":"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.","section":"Section 3, coverage assessment"}],"minor_comments":[{"comment":"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.","section":"Section 2.3.2"},{"comment":"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.","section":"Algorithm 1, line 10"},{"comment":"The caption contains a typo: 'the Gaussian mixture model g(epsilon) give in (4)' should read 'given in (4)'.","section":"Figure 1 caption"},{"comment":"The phrase 'In such a case case' contains a duplicated word and should be 'In such a case'.","section":"Remark 4"},{"comment":"The GIG parameter in the update for lambda is typeset as '-n + 2/2', which is ambiguous; it should be '-(n-2)/2'.","section":"Section 2.3.1"},{"comment":"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.","section":"Section 3"}],"recommendation":"major_revision","confidential_remarks":"The paper's main contribution is a Bayesian algorithmic framework for discretization error quantification. The framework is potentially useful, but the printed Gibbs sampler does not target the intended posterior because of the incorrect s_i update, and the tau_1/tau_j typo compounds the problem. These are fixable, but the numerical experiments must be rerun with a corrected sampler, and the empirical validation should be strengthened with repeated-simulation coverage rather than single-trajectory in-sample plots. I recommend major revision rather than rejection, as the conceptual approach is sound and the errors are local. There are no citation or scope concerns; the fit with the journal is appropriate."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"The stress-test note is right, and it's the first thing you should know: the Gibbs sampler in Section 2.3.2 and Algorithm 1 does not target the stated posterior. The latent mixture indicators s_i are drawn from the marginal prior (5), but a correct Gibbs step must sample them from the conditional posterior given the current η and the data, p(s_i = k | z_i, η_i) ∝ w_k N(z_i − η_i; m_k, v_k^2). Sampling from the prior breaks the data augmentation. The numerical sections therefore don't demonstrate the method as advertised. There's also a typo in the variance update for η_j: the formula w_j^2 = \tilde{w}_j^2 λτ_1/(\tilde{w}_j^2 + λτ_j) should have τ_j in the numerator, not τ_1, since the prior precision is 1/(λτ_j).\n\nWhat's actually new and good: extending Bayesian isotonic regression from normal means to normal variances, using the log-χ² Gaussian mixture and a horseshoe-type shrinkage prior on log increments. That's a natural idea for probabilistic numerics, and the conditional derivations for η_j given s are correct in spirit. If the sampler were fixed, the method would likely work and fill a specific niche.\n\nOther soft spots are secondary but real: the validation is in-sample, with credible intervals computed from the same residuals used to fit; no code is released for the new method; and the assumptions of independence, Gaussianity, and monotone variance growth are acknowledged in Remark 2 but not stress-tested. The paper also doesn't compare with existing ODE filters or probabilistic integrators.\n\nThis paper is for people in probabilistic numerics or inverse problems who want a Bayesian handle on discretization-error variance. It deserves a serious referee, not a desk reject, but it needs substantive revision. I'd condition acceptance on a corrected sampler and released code, and I'd want to see the numerical results rerun before believing the coverage claims.","headline":"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.","tokens_in":13150,"tokens_out":3058,"would_cite":false,"duration_ms":32357,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["65L05","62F15","62G08"],"pacs":[],"model":"deepseek-v4-flash","headline":"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…","keywords":["discretization error quantification","Bayesian isotonic regression","horseshoe prior","Gibbs sampler","Gaussian mixture model","ordinary differential equations","probabilistic numerics","uncertainty quantification"],"falsifier":"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.","tokens_in":12163,"feed_emoji":"📊","tokens_out":7405,"duration_ms":64331,"temperature":0.7,"pith_summary":"This paper develops a Bayesian way to put error bars on the numerical solution of an ordinary differential equation when only noisy observations of the true solution are available. The key idea is to treat the discretization error at each observation time as a zero-mean Gaussian whose variance grows monotonically as integration proceeds, and to estimate those variances from the residuals between observations and numerical approximations. The authors show that a horseshoe-type shrinkage prior on the log-increments of the variances, together with a 10-component Gaussian mixture approximation of the log-chi-square distribution, turns the posterior into a Gibbs sampler whose full conditionals are all standard distributions. In numerical tests on the FitzHugh–Nagumo and Kepler systems, the resulting credible intervals for the absolute residual and error capture the observed scale of the error. If the framework is right, it gives practitioners a principled way to quantify and propagate discretization uncertainty in Bayesian inverse problems for ODEs.","feed_headline":"Posterior intervals track ODE discretization error","feed_subtitle":"A horseshoe prior and 10-component Gaussian mixture turn error-variance estimation into a standard Gibbs sampler.","key_machinery":"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.","core_discovery":"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.","pith_inferences":["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."],"forward_implications":["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."],"supporting_citations":[{"why":"Supplies the half-shrinkage prior construction and hierarchical representation that the paper adapts from normal means to normal variances.","marker":"[16]"},{"why":"The previous maximum-likelihood isotonic approach for discretization error quantification that motivates the model and serves as the comparison baseline in the experiments.","marker":"[13]"},{"why":"Provides the 10-component Gaussian mixture weights, means, and variances used to approximate the log-chi-square distribution.","marker":"[29]"},{"why":"Introduces the mixture-sampler and data-augmentation idea that the paper's Gibbs sampler builds on.","marker":"[28]"},{"why":"Earlier work modelling discretization error with the Wishart distribution, showing the line of work and the context for the diagonal-variance restriction.","marker":"[12]"},{"why":"Origin of the log-variance transformation used to turn variance parameters into location parameters in a regression-like likelihood.","marker":"[27]"},{"why":"Provides the generalized inverse Gaussian random-variate algorithm needed for stable sampling of the full conditionals in the Gibbs sampler.","marker":"[25]"}],"fun_headline_variants":["Bayesian isotonic regression quantifies ODE solver error","Horseshoe prior tracks discretization error in ODE solvers","Bayesian monotone model quantifies ODE integration error","Posterior tracks ODE solver error via isotonic regression","Gibbs sampler enables Bayesian error estimation for ODE solvers"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"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.","fun_headline_variants_meta":{"raw":{"variants":["Bayesian isotonic regression quantifies ODE solver error","Horseshoe prior tracks discretization error in ODE solvers","Bayesian monotone model quantifies ODE integration error","Posterior tracks ODE solver error via isotonic regression","Gibbs sampler enables Bayesian error estimation for ODE solvers"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000867,"raw_usage":{"total_tokens":3769,"prompt_tokens":970,"completion_tokens":2799,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":586,"completion_tokens_details":{"reasoning_tokens":2729}},"tokens_in":586,"tokens_out":2799,"duration_ms":18056,"temperature":1.0,"reasoning_tokens":2729,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-12T21:42:13.237132+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"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.","supporting_citations":[{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"The previous maximum-likelihood isotonic approach for discretization error quantification that motivates the model and serves as the comparison baseline in the experiments."},{"cited_title":"PhD thesis, Massachusetts Institute of Technology (1988) 18","cited_arxiv_id":null,"evidence_quote":"Origin of the log-variance transformation used to turn variance parameters into location parameters in a regression-like likelihood."},{"cited_title":"Diploma Thesis, Technical University Graz, Austria (1989)","cited_arxiv_id":null,"evidence_quote":"Provides the generalized inverse Gaussian random-variate algorithm needed for stable sampling of the full conditionals in the Gibbs sampler."}],"review_version":1}