Pith. sign in

REVIEW 3 major objections 4 minor 22 references

For weight functions outside classical families, statisticians can compute custom Gauss quadrature rules to about 10^-18 (moment determinants) or 10^-15 (Stieltjes) accuracy using two elementary methods, with an R-callable implementation.

Reviewed by Pith at T0; open to challenge. T0 means a machine referee read the full paper against a public rubric. the ladder, T0–T4 →

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗

CustomGaussQuadrature computes custom Gauss quadrature nodes and weights for non-classical weight functions to roughly 10^-18 (moment determinants) or 10^-15 (Stieltjes) accuracy, using adaptive-precision arithmetic and an R interface.

T0 review reviewed 2026-08-02 challenge →

load-bearing objection Honest software paper for custom Gauss rules, but the accuracy certificates are based on internal self-agreement and can certify a wrong Stieltjes rule on narrow peaks. the 3 major comments →

arxiv 2607.14511 v1 pith:RY7OWP4U submitted 2026-07-16 stat.CO

Custom-made Gauss quadrature: an introduction for statisticians

Paul Kabaila This is my paper
classification stat.CO MSC 65D3265D3033C4565-04
keywords Gauss quadraturecustom-made Gauss rulesmoment determinantsStieltjes procedureorthogonal polynomialshigh-precision arithmeticJuliaR interface
verification ladder T0 review T1 audit T2 compute T3 formal T4 reserved

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

Gauss quadrature rules are exact for polynomials of degree at most 2n-1 and need only n function evaluations, but ready-made rules exist only for classical weights such as Legendre, Hermite, or Laguerre. This paper argues that for any other weight function, a statistician can still obtain a near-double-precision rule by computing the recursion coefficients of the associated monic orthogonal polynomials, either from moments via Hankel determinants or from the weight via the Stieltjes procedure. The paper contributes an implementation that adaptively raises the arithmetic precision and, for the Stieltjes route, the size of an auxiliary Gauss-Legendre rule, until two successive computations agree, and it treats that agreement as a certified error bound. If correct, statisticians can accurately compute expected values under non-classical densities such as scaled chi, Weibull, and inverse gamma without leaving R. The motivating applications include simultaneous inference and multivariate t probabilities, where the integrand is itself expensive to evaluate.

Core claim

The central claim is that custom-made Gauss rules can be computed to about double-precision accuracy without specialized knowledge. Step 1 computes the coefficients alpha_k and beta_k in the three-term recurrence for monic orthogonal polynomials, either from the weight's moments using Hankel-determinant formulas or from the weight itself using the Stieltjes procedure with an auxiliary Gauss-Legendre rule. Step 2 turns those coefficients into nodes and weights by diagonalizing the symmetric tridiagonal Jacobi matrix. The novel part is the error control: BigFloat arithmetic at progressively higher bit counts (up to 444 bits for moment determinants; 256-512 bits for the Stieltjes route) is run

What carries the argument

The carrying mechanism is the three-term recurrence pi_{k+1}(x) = (x - alpha_k) pi_k(x) - beta_k pi_{k-1}(x), with coefficients alpha_k = (x pi_k, pi_k) / (pi_k, pi_k) and beta_k = (pi_k, pi_k) / (pi_{k-1}, pi_{k-1}), which determines the monic orthogonal polynomials and hence the Gauss rule. Two routes feed it: moment determinants, which express alpha and beta as ratios of Hankel determinants in the moments mu_s, and the Stieltjes procedure, which builds the same coefficients iteratively by approximating inner products with an r-node auxiliary Gauss-Legendre rule after mapping the support interval to (-1,1). The error control is carried by the self-agreement rule: compute the same ill-condi

Load-bearing premise

The error certificate is the agreement between two finite-precision computations of the same map; if both precisions are too low together (e.g., moments overflow even 444-bit arithmetic, or the auxiliary Stieltjes rule converges smoothly but slowly), the chosen difference can be tiny while the rule is far from accurate.

What would settle it

Use the package to compute the n=6 Gauss rule for the Weibull weight with k=3.1, then integrate g(x)=x^12 (a degree within the exactness range) and compare with Gamma(1 + 12/3.1) evaluated in 1000-digit arithmetic. A relative error above about 10^-15 would show the self-agreement certificate is not a valid accuracy bound. Equally, feed a classical weight (e.g., the Legendre weight on (-1,1)) through the non-classical path; the returned rule must match the known Legendre nodes and weights to the stated tolerance.

Watch this falsifier. Get emailed when new claim-graph text bears on it.

If this is right

  • For any weight function for which a moment formula is supplied, the moment-determinants route yields a Gauss rule with assessed maximum absolute node error and maximum relative weight error below 10^-18.
  • For any weight function for which log f and mu0 are supplied, the Stieltjes route yields the same quantities below 10^-15, at the cost of adaptively choosing the auxiliary rule size r, which can be surprisingly large (e.g., Weibull with small shape parameter).
  • Once the recursion coefficients are known, rules with 1,2,...,n nodes follow from the same coefficients, so one computation serves many quadrature sizes.
  • An R user can obtain these rules without leaving R, using the Julia interface described in the paper, by supplying a short snippet for the moments or for log f.
  • For distributions on unbounded supports, the initial transformation of the support interval to (-1,1) and the log-domain computation of quadrature weights avoid overflow and NaN issues that would otherwise arise.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • If the internal self-agreement certificate is valid, the same adaptive-precision pattern could be applied to other ill-conditioned moment-based computations, such as the construction of quadrature rules with prescribed nodes or the computation of moments of truncated distributions; the paper itself does not test this transfer.
  • The stated 10^-18 and 10^-15 bounds are not cross-checked against an independent external reference; a user wanting certifiable accuracy should verify against a known integral before trusting the bounds in a safety-critical setting.
  • The moment-determinants route has a hard ceiling of 444 bits of precision; weight functions whose moments fail to be stably computed within that range may fall outside the package's guaranteed regime, even though the Stieltjes route may still work.
  • For the Stieltjes route, the auxiliary rule size is chosen by successive r values; a slowly converging discretization could in principle satisfy the agreement test at a moderate r while still being far from the true rule, so a pragmatic check is to also compare against the moment-determinants rule when both routes are available.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit.

Referee Report

3 major / 4 minor

Summary. The paper describes and implements two classical methods for computing custom n-point Gauss quadrature rules for non-classical weight functions: the moment-determinant method and the Stieltjes (discretization) procedure. The main practical contribution is the Julia package CustomGaussQuadrature, which uses type-generic programming, adaptive BigFloat precision, and a self-agreement rule to certify that the computed recursion coefficients have errors bounded by about 10^-18 (moment determinants) or 10^-15 (Stieltjes). Access from R via JuliaConnectoR is illustrated on Weibull and inverse-gamma weight functions. The underlying mathematics is standard, but the paper's accuracy certification is internal and the numerical evidence consists only of code snippets rather than reported results.

Significance. If the stated accuracy guarantees are reliable, the package would be a useful, easy-to-call tool for statisticians needing custom Gauss rules, complementing existing packages such as ORTHPOL and PolyChaos. The paper is honest about limitations: Appendix A states there is no guarantee that the moment-determinant bootstrap succeeds, and Appendix B allows rejection when precision agreement fails. The two methods cross-check each other on the two examples, which is a positive feature. However, the central claim is not yet independently verified: the accuracy certificate is based on comparing two computations from the same package, and the Stieltjes r-selection test can certify a wrong rule when the transformed integrand has an unresolved narrow peak. These issues are fixable within the scope of the manuscript, so the paper warrants a major revision rather than rejection.

major comments (3)
  1. [§4.1 and Appendix A] The moment-determinant accuracy certificate is entirely internal. Equations (11)–(12) accept the computation when two BigFloat precisions b_j and b_{j+1} agree to 1e-18, and the Introduction states that |c1 - c2| is treated as an upper error bound. This presupposes that the higher-precision computation is much more accurate than the lower. Nothing in the paper rules out a systematic error that is nearly identical at both precisions, e.g., a moment formula that is inaccurate for the chosen parameter range, or determinants that saturate the cap b ≤ 444 in a way that cancels in the difference. The method is also validated only against the same package's Stieltjes implementation, so the two-method agreement is not an independent check. Please validate against an external standard (e.g., Gautschi's ORTHPOL) or against integrals with known exact values, and report the actual differences.
  2. [Appendix B (choice of r)] The r-selection rule can stop prematurely on an unresolved narrow peak. Candidate sizes are r_j = j(n+7), and r_final is the first r_{j+1} for which the recursion coefficients at r_j and r_{j+1} agree to 1e-15. If h(y) in Eq. (9) has a peak narrower than the Gauss–Legendre spacing at both r_j and r_{j+1}, both discrete rules can miss the peak and return nearly identical wrong inner products; the difference test then passes and the algorithm certifies an inaccurate rule. The paper's own Weibull example (Sections 4.2 and 6.3) has exactly this behaviour as k → 0, with the integrand peak width approaching 0. The paper notes that r increases as k decreases, but it does not acknowledge that the r-difference test is not a valid certificate during the unresolved phase. Please add a safeguard, for example requiring agreement with the moment-determinant rule for small k, checking moment identities
  3. [§6 (Examples 1–4)] The manuscript contains no reported numerical results. It shows code snippets and instructs the reader to run cat(...) comparisons between the moment-determinant and Stieltjes outputs, but it never reports the resulting nodes, weights, maximum absolute differences, or timings. As a result, the central quantitative claims — 'maximum absolute errors ... bounded above by 10^-18' (Section 4.1) and 'about 10^-15' (Section 4.2) — are not demonstrated in the paper itself. A compact table for the Weibull and inverse-gamma examples, including a check against a known exact integral, would make the claims assessable and would substantially strengthen the paper.
minor comments (4)
  1. [§4.2] Typo: 'we we use' should read 'we use'. Also, 'for a ga polynomial of degree 2n-1' should be 'for a polynomial of degree 2n-1'.
  2. [§3] In the formula for the weight λ_i, the symbol x_i is used both for the node and for the eigenvector, which is confusing. Use a distinct notation, e.g., v_i for the eigenvector.
  3. [Appendix B] The text first says 'All that the user needs to provide is a Julia function to evaluate f', but then explains that the user must provide log(f) to avoid NaN issues. Please make the final requirement explicit at the beginning of the appendix.
  4. [References] Some references, such as the PolyChaos and JuliaConnectoR versions, are dated 2026; ensure the list is consistent and final. Also, use the accent in 'Fejér' consistently.

Circularity Check

0 steps flagged

No significant circularity: the quadrature computation is grounded in external standard theory; internal self-agreement checks are fallible validation heuristics, not definitional reductions.

full rationale

The paper's derivation chain is standard and externally grounded: user-supplied moments or log-weight are used to compute recursion coefficients (moment determinants per Gautschi 2004 Theorem 2.2, or the Stieltjes discretization per Gautschi 2004 §§2.2.2–2.2.3), and the nodes/weights are then obtained from the Jacobi matrix eigensystem (Section 3, with conditioning results from Stewart 2001). None of these steps defines the nodes/weights in terms of the accuracy checks. The adaptive-precision and r-selection rules (Appendix A stopping when (11) and (12) are below 10^-18; Appendix B choosing r_final when successive a_i and b_i differ by at most 10^-15) are self-referential quality checks: they assert that two high-precision computations agree and then take that agreement as a practical error estimate. This is explicitly acknowledged by the paper: 'we take c2 as our final approximation and regard |c1 − c2| as a practical upper estimate of its absolute error' (Introduction). That is a fallible heuristic, not a circular derivation: the final c2 is not constructed from c1, and the 10^-18 / 10^-15 statements are presented as assessments ('expected to have', 'assessed by this method') rather than as theorems. The same applies to the Stieltjes r-difference test: the narrow-peak scenario for Weibull k→0 is a genuine convergence risk, and the paper itself flags that 'surprisingly large values of r may need to be chosen' (Section 6.3), but that is a correctness risk, not circularity, because the stopping rule does not define the quadrature rule. The self-citations (Kabaila 2021, 2022, 2023, 2026) are contextual, historical, or package references; they are not load-bearing evidence for the recurrence/eigenvalue mathematics, which is drawn from Gautschi and Stewart. I therefore find no step that reduces to its own input by construction.

Axiom & Free-Parameter Ledger

5 free parameters · 6 axioms · 0 invented entities

The central claim rests on textbook mathematics (Gautschi's theorems, Stewart's eigensystem results) plus five hand-chosen algorithmic constants — precision ladders, agreement thresholds, and the auxiliary-rule grid — that determine the claimed accuracy. The load-bearing domain assumptions are that Stieltjes discretization converges regularly in r and that two-precision self-agreement bounds the true error; both are asserted and motivated but not externally validated. No invented entities: the package is a software artifact, not a postulated entity. Note that the headline accuracy numbers (10^-18, 10^-15) are exactly the agreement thresholds, so the reported accuracy is an input target of the certification scheme rather than an independently measured output.

free parameters (5)
  • moment-determinants precision stepping = 26-bit increments, max 444 bits
    Appendix A: b steps through b1, b1+26, ... up to 444; the error estimate compares successive steps. Hand-chosen ladder.
  • moment-determinants accuracy target = 10^-18
    Appendix A: b_final is the first precision where successive α and β differences drop below 10^-18; this is the claimed accuracy, set as an input target rather than measured independently.
  • Stieltjes auxiliary-rule grid = r_j = j(n+7), j = 3..40 (defaults)
    Appendix B: offset 7 and j_max 40 bound the auxiliary Gauss-Legendre rule sizes tested; if convergence has not appeared by j_max, the problem is rejected rather than solved.
  • Stieltjes agreement thresholds = 10^-22 (precision agreement), 10^-15 (r convergence)
    Appendix B: these thresholds determine when the procedure stops and what accuracy is reported.
  • Stieltjes precision ladder = 128/224/256/480/512 bits
    Appendix B: fixed candidate working precisions compared pairwise; 128-bit emulated by Double64. The final certified accuracy is a byproduct of this hand-chosen set.
axioms (6)
  • standard math The moment condition ∫|x|^{2n−1} f dx < ∞ guarantees existence of the n-point Gauss rule
    Section 1, stated as the existence condition for the rule.
  • standard math Monic orthogonal polynomials satisfy the three-term recurrence (3) with α_k, β_k from (4) (Gautschi 2004, Thm 1.27)
    Section 2.1, the basis of both Step-1 methods.
  • standard math Recursion coefficients are given by the Hankel-determinant ratios (5)-(6) (Gautschi 2004, Thm 2.2)
    Section 2.2, the moment-determinant method.
  • domain assumption The Stieltjes procedure's recursion coefficients converge to the exact values as the auxiliary-rule size r → ∞
    Section 2.3 says convergence is 'expected'; Appendix B's stopping rule assumes the convergence is regular enough for the 10^-15 self-agreement test to be reliable.
  • domain assumption Agreement between two finite-precision computations (|c1 − c2|) is a practical upper bound on the true error
    Introduction and Appendices A/B; load-bearing for every reported error bound and never checked against an external standard.
  • standard math Eigenvalues of the Hermitian Jacobi matrix are perfectly conditioned, and eigenvectors are well-conditioned when nodes are well separated (Stewart 2001)
    Section 5; justifies that Step 2 preserves the accuracy of the recursion coefficients.

reviewed 2026-08-02 · how reviews work

0 comments
Cite this review

Pith. "Pith review of Custom-made Gauss quadrature: an introduction for statisticians." pith.science (2026). https://pith.science/paper/RY7OWP4U

@misc{pith2026260714511,
  author       = {Pith},
  title        = {Pith review of: Custom-made Gauss quadrature: an introduction for statisticians},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/RY7OWP4U}},
  note         = {Machine review of arXiv:2607.14511}
}
Share X Bluesky LinkedIn Reddit HN
abstract

An $n$-point Gauss quadrature rule approximates the weighted integral of a function by a weighted average of $n$ evaluations of this function and is exact for polynomials of degree at most $2n-1$. Such rules can be highly accurate with relatively few evaluations. For weight functions that are associated with classical orthogonal polynomials of a continuous variable (such as Legendre, Hermite and Laguerre), these rules are readily available. We suppose that this is not the case, so that these rules must be custom-made. The two most easily understood methods for the computation of these rules are (a) moment determinants and (b) the Stieltjes procedure. We implement them in the Julia package CustomGaussQuadrature, which uses type-generic numerical programming and adaptive high-precision arithmetic to assess the approximation error due to roundoff. We describe access from R via JuliaConnectoR.

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

22 extracted references

  1. [1]

    Chihara, T.S. 1978. An Introduction to Orthogonal Polynomials . Mathematics and Its Applications, vol. 13. New York: Gordon and Breach

    1978

  2. [2]

    Dunnett, C.W. 1989. Algorithm AS 251: Multivariate normal probability integrals with product correlation structure. Journal of the Royal Statistical Society, Series C (Applied Statistics) , 38 , 564--579

    1989

  3. [3]

    Dunnett, C.W., & Sobel, M. 1955. Approximations to the probability integral and certain percentage points of a multivariate analogue of S tudent's t-distribution. Biometrika , 42 , 258--260

    1955

  4. [4]

    Gautschi, W. 1968. Construction of G auss- C hristoffel quadrature formulas. Mathematics of Computation , 22 , 251--270

    1968

  5. [5]

    Gautschi, W. 1983. How and how not to check G aussian quadrature formulae. BIT , 23 , 209--216

    1983

  6. [6]

    Gautschi, W. 1994. Algorithm 726: ORTHPOL - A package of routines for generating orthogonal polynomials and G auss-type quadrature rules. ACM Transactions on Mathematical Software , 20 , 21--62

    1994

  7. [7]

    Gautschi, W. 2004. Orthogonal Polynomials: Computation and Approximation . New York: Oxford University Press

    2004

  8. [8]

    Genz, A., & Bretz, F. 2009. Computation of Multivariate Normal and t Probabilities . London: Springer

    2009

  9. [9]

    Gupta, S.S., & Panchapakesan, S. 2002. Multiple Decision Procedures: Theory and Methodology of Selecting and Ranking Populations . Philadelphia: SIAM

    2002

  10. [10]

    Hochberg, Y., & Tamhane, A.C. 1987. Multiple Comparison Procedures . New York: Wiley

    1987

  11. [11]

    Johnson, N.L., Kotz, S., & Balakrishnan, N. 1995. Continuous Univariate Distributions . second edn. Vol. 2. New York: John Wiley

    1995

  12. [12]

    Johnson, Steven G. 2026. QuadGK.jl : G auss-- K ronrod integration in J ulia . https://github.com/JuliaMath/QuadGK.jl. V ersion 2.11.3

    2026

  13. [13]

    Kabaila, P. 2022. custom.gauss.quad : R package for custom-made G auss quadrature . https://CRAN.R-project.org/package=custom.gauss.quad. V ersion 1.0.0

    2022

  14. [14]

    Kabaila, P. 2026. CustomGaussQuadrature.jl : Julia package for custom-made G auss quadrature . https://github.com/pvkabaila/CustomGaussQuadrature.jl. V ersion 3.0.6

    2026

  15. [15]

    Kabaila, P., & Perera, A. 2023. Model averaged tail area confidence intervals in nested linear regression models. Australian & New Zealand Journal of Statistics , 65 , 364--378

    2023

  16. [16]

    Kabaila, P., & Ranathunga, N. 2021. Computation of the expected value of a chi-distributed random variable. Computational Statistics , 36 , 313--332

    2021

  17. [17]

    Lenz, S., Binder, H., D'Ambrosio, A., & Choe, J. 2026. Julia C onnecto R : A Functionally Oriented Interface for Integrating ' J ulia' with R . https://CRAN.R-project.org/package=JuliaConnectoR. V ersion 1.1.5

    2026

  18. [18]

    Miller, R.G. 1981. Simultaneous Statistical Inference . 2nd edn. New York: Springer

    1981

  19. [19]

    Miwa, T., Hayter, A.J., & Kuriki, S. 2003. The evaluation of general non-centred orthonant probabilities. Journal of the Royal Statistical Society, Series B , 65 , 223--234

    2003

  20. [20]

    Mori, M. 1988. The double exponential formula for numerical integration over the half infinite interval. Pages 367--379 of: Agarwal, R.P., Chow, Y.M., & Wilson, S.J. (eds), Numerical Mathematics (Singapore 1988) . Basel: Birkhauser

    1988

  21. [21]

    M \"u hlpfordt , Tillmann, Zahn , Frederik, Hagenmeyer , Veit, & Faulwasser , Timm. 2026. PolyChaos.jl : Orthogonal Polynomials, Quadrature, and Polynomial Chaos . https://github.com/SciML/PolyChaos.jl. V ersion 1.2.0

    2026

  22. [22]

    Stewart, G.W. 2001. Matrix Algorithms. Volume II: Eigensystems . Philadelphia: Society for Industrial and Applied Mathematics

    2001

This paper was first reviewed by deepseek-v4-flash on August 2, 2026.