{"id":"99447e05-9ab2-416f-93bd-4ac8962af8db","arxiv_id":"2607.14511","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":5.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":5,"one_line_summary":"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.","lead":"Paul Kabaila's paper presents a Julia package, CustomGaussQuadrature, that computes Gauss quadrature rules for non-classical weight functions using two classic methods armored with adaptive high-precision arithmetic. Statisticians who need accurate integrals under non-standard densities such as Weibull or inverse gamma can call the package from R and obtain rules whose error is certified by self-agreement checks.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Stieltjes r-selection can stop prematurely on a narrow unresolved peak, so the 1e-15 self-agreement certificate may certify a wrong rule (e.g., Weibull k→0).","rationale":"The reader's weakest assumption was that the self-agreement error test is a valid certificate of accuracy. My concern sharpens that: the Stieltjes r-selection rule has a concrete failure mode in which agreement between successive r values is achieved by both rules under-resolving the same narrow feature. This is load-bearing because the central 1e-15 accuracy claim for the Stieltjes path rests on that test. I do not change the verdict because the paper is candid about the difficulty (Section 6.3), and the moment-determinant path is not affected. The fix is internal: before stopping, verify polynomial exactness against the known moments (or compare with the moment method); the proposed test checks whether the current implementation already does this. If the test passes for small k, the concern is laid to rest. If it fails, the paper's stated error bounds need qualification. The concern is not about external validation per se, but about an internal flaw in the stopping rule that external validation would uncover.","tokens_in":13467,"tokens_out":10990,"duration_ms":112548,"concrete_test":"Run the package's Stieltjes procedure for the Weibull pdf with k=0.1, n=6, default settings. If it returns a rule, compute max_{m=0..11} |Σ_i λ_i τ_i^m − Γ(1+m/k)| / Γ(1+m/k). If this exceeds 1e-10, the self-agreement certificate has certified a rule whose error is far above the claimed 1e-15. Repeat for k=0.05 and k=0.01, and record whether the procedure rejects the problem or silently returns an inaccurate rule.","verdict_should_be":"UNCHANGED","load_bearing_attack":"In Appendix B, the Stieltjes procedure selects r_final = r_{j+1} for the smallest j such that the recursion coefficients computed at r_j and r_{j+1} agree to 1e-15 (relative for sqrt(beta)). The candidate r_j = j(n+offset) increases by only offset=7 each step. If the transformed integrand h(y) in (9) has a very narrow peak whose width is smaller than the Gauss-Legendre node spacing at both r_j and r_{j+1}, both discrete rules miss the peak and return nearly identical (wrong) inner products. The difference test then passes, and the algorithm stops with an inaccurate rule without ever resolving the feature. This is not a hypothetical: for the paper's own Weibull example, the integrand for moment s has a peak whose width shrinks rapidly as k→0 (Section 4.2 and 6.3). The paper acknowledges that r increases as k→0, but does not state that the r-difference test can certify a wrong rule during the unresolved phase. The moment-determinant method is unaffected because moments are evaluated analytically; the discrepancy would be visible if the two methods were compared, but the paper presents no such comparison for small k.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","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.","tokens_in":13642,"tokens_out":4486,"duration_ms":48277,"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":[{"comment":"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.","section":"§4.1 and Appendix A"},{"comment":"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","section":"Appendix B (choice of r)"},{"comment":"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.","section":"§6 (Examples 1–4)"}],"minor_comments":[{"comment":"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'.","section":"§4.2"},{"comment":"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.","section":"§3"},{"comment":"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.","section":"Appendix B"},{"comment":"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.","section":"References"}],"recommendation":"major_revision","confidential_remarks":null},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Here's my take. The paper is an honest, well-written software paper. It restates two standard methods from Gautschi (2004) — moment determinants and the Stieltjes procedure — for computing Gauss quadrature for non-classical weight functions, and contributes a Julia package CustomGaussQuadrature with adaptive high-precision arithmetic and an R interface via JuliaConnectoR. The math is standard and correctly described; the code is available and the examples run. The paper is candid about limitations: it notes that the Stieltjes procedure can require very large r, that the moment-method bootstrap 'has no guarantee' (Appendix A), and that the error estimate is a 'practical upper estimate' rather than a rigorous bound. That honesty is genuine and should be credited.\n\nThe soft spot is the accuracy certification. The advertised 10^-18 and 10^-15 error bounds are based on agreement between two computations at different precisions (or different auxiliary rule sizes). This is an internal consistency check, not an external validation. It can, in principle, certify a wrong answer. The stress-test note makes a specific and plausible case for the Stieltjes procedure: r increases in steps of offset=7, and the algorithm stops when two consecutive r values agree to 1e-15. If the transformed integrand h(y) has a peak narrower than the Gauss-Legendre node spacing at both r values, both rules miss the peak and agree on the same wrong inner product. The paper's own Weibull example has this kind of peak-shrinking behavior as k→0; the paper acknowledges that r grows but does not mention the false-convergence failure mode. So a user with a sharply-peaked weight function could get a certified-looking but inaccurate rule. The moment-determinant method is safer because the moments are evaluated analytically, but the package offers no external cross-check against, say, ORTHPOL or an integral with a known exact value, and only n=6 and n=8 examples are shown.\n\nIn proportion: the moment-method implementation seems solid for the examples, and the Stieltjes method's risk is a real but specific failure mode, not a general collapse. The paper's central pedagogical value holds up. What's missing is a demonstration that the self-agreement test actually predicts accuracy on a problem with a known answer. That is easy to add.\n\nWho is this for? Statisticians who need custom Gauss rules and want a friendly, R-callable implementation, plus readers who want a compact introduction to the two methods. It deserves serious refereeing. My recommendation: send it to review, and have the authors add an ORTHPOL comparison, an exact-integral test, and a large-n or extreme-parameter case. With that, the accuracy claims would be on much firmer ground.","headline":"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.","tokens_in":14284,"tokens_out":4898,"would_cite":false,"duration_ms":43638,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["65D32","65D30","33C45","65-04"],"pacs":[],"model":"deepseek-v4-flash","headline":"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.","keywords":["Gauss quadrature","custom-made Gauss rules","moment determinants","Stieltjes procedure","orthogonal polynomials","high-precision arithmetic","Julia","R interface"],"falsifier":"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.","tokens_in":13206,"feed_emoji":"🧮","tokens_out":9136,"duration_ms":78429,"temperature":0.7,"pith_summary":"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.","feed_headline":"Custom Gauss rules to 18-digit accuracy from an R call","feed_subtitle":"Adaptive high-precision arithmetic certifies nodes and weights to 10^-15–10^-18.","key_machinery":"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","core_discovery":"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","pith_inferences":["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."],"forward_implications":["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."],"fun_headline_variants":["Custom Gauss quadrature, certified to 18 digits","Precision-controlled Gauss rules for any weight","From R: tailor-made Gauss rules, error-certified","Adaptive precision Gauss quadrature, R-friendly","Custom Gauss nodes and weights, accuracy assured"],"cache_read_input_tokens":2304,"weakest_assumption_plain":"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.","fun_headline_variants_meta":{"raw":{"variants":["Custom Gauss quadrature, certified to 18 digits","Precision-controlled Gauss rules for any weight","From R: tailor-made Gauss rules, error-certified","Adaptive precision Gauss quadrature, R-friendly","Custom Gauss nodes and weights, accuracy assured"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000147,"raw_usage":{"total_tokens":990,"prompt_tokens":676,"completion_tokens":314,"prompt_tokens_details":{"cached_tokens":256},"prompt_cache_hit_tokens":256,"prompt_cache_miss_tokens":420,"completion_tokens_details":{"reasoning_tokens":242}},"tokens_in":420,"tokens_out":314,"duration_ms":3292,"temperature":1.0,"reasoning_tokens":242,"cache_read_input_tokens":256,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-02T01:52:41.088624+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"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.","supporting_citations":[],"review_version":1}