{"id":"ea00bc45-91e2-4c33-a4f7-84f50279f79d","arxiv_id":"2411.12938","paper_version":2,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":5.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":0,"one_line_summary":"The authors benchmark and release two numerical routes to ratio-distribution PDFs, a 1D double-exponential Mellin convolution and a 2D vectorized Broda-Khan characteristic-function inversion, applied to Hake's ratio and one non-normal case.","lead":"This paper benchmarks two numerical ways to compute the probability distribution of a ratio of random variables, using Hake's gain in education research as the test case. It matters because exact ratio distributions are needed in metrology and statistics, where Monte Carlo is slow and closed-form formulas often do not exist.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Theorem 6's Broda-Kan PDF integral is not absolutely convergent for the tested normal cases, and Appendix C's algorithm appears to omit the second term of the difference form, so the 2D method's reported accuracy lacks a valid theoretical basis.","rationale":"The central claim covers two computational methods. The 1D DE-Mellin approach is well supported: it is benchmarked against PaCAL and against analytic reference values computed in quadruple precision, so I have no significant objection to that part. The 2D Broda-Khan inversion is the load-bearing weak point. The reader's weakest assumption flagged the unverified absolute-convergence condition in Theorem 6. My reading sharpens this: for the normal-normal case, the integrand in (B4) has a non-removable 1/t singularity, so the integral cannot converge absolutely for the tested parameters; the theorem's hypothesis is not merely unverified but false as stated. Moreover, the algorithm in Appendix C integrates t2 over positive values only, with a single term, which does not match either the principal-value form of (B4) or the difference form obtained by differentiating the CDF. This is an internal inconsistency in the manuscript. If the actual code implements the difference form or includes a factor of 2, then the paper's formulas and pseudo-code need correction; if the code follows the appendix literally, the reported 2D accuracy is not trustworthy for the pilot cases. The standard-normal ratio test at x=0 would settle which situation holds. Because the paper's overall contribution also includes the solid 1D method and the statistical review of Hake's ratio, this concern does not overturn the paper entirely, but it does reinforce the need for the CONDITIONAL verdict already given. I therefore recommend UNCHANGED: the reader's conditional acceptance remains appropriate, with the added requirement that the authors clarify Theorem 6 and verify the 2D implementation on a closed-form test case.","tokens_in":20433,"tokens_out":18355,"duration_ms":173422,"concrete_test":"Run the published GitHub code, or re-implement Algorithm 1 exactly as written (single sum over t2>0, c=h1h2/π², integrand φ_2(-(x t1+t2))φ_1(t1)/t2), for the standard normal ratio case (a,b)=(0,0) at x=0; the exact Cauchy PDF value is 1/π, so an output of 1/(2π) would confirm that the discretization omits the second term and does not compute the PDF of (3.13).","verdict_should_be":"UNCHANGED","load_bearing_attack":"Appendix B's Theorem 6 states that the 2D PDF inversion formula (3.13)/(B4), f_R(r)=π^{-2}∫_0^∞∫_{-∞}^{∞} Re[φ_2(s,-t-rs)]/t ds dt, holds 'whenever this integral converges absolutely.' For the pilot's normal-normal case, φ_2(s,τ)=φ_{X1}(s)(ib-τ)e^{ibτ-τ²/2}; at τ=-t-rs, Re[φ_2]/t has a non-removable 1/t singularity at t=0 with coefficient rs·cos(brs)+b·sin(brs), which is nonzero for the tested (a,b) pairs. Hence the integral is not absolutely convergent and the theorem's hypothesis is not satisfied. Differentiating the CDF (B3) instead yields the difference form over t>0, which has a removable singularity and is a different integral. Appendix C's Algorithm 1 sums only t2>0 of the single term φ_2(-(x t1+t2)), omitting the φ_2(s,t-xs) contribution; for (a,b)=(0,0) this would give half the Cauchy density at x=0. The paper does not state a principal-value interpretation, verify absolute convergence, or reconcile the theorem with the algorithm. Therefore the 2D method's reported accuracy is not linked to a valid stated theorem, and the claim of a reliable, reproducible 2D implementation is insecure without a corrected statement and an implementation check.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper proposes two numerical methods for computing the density of a ratio of random variables: a one-dimensional double-exponential (DE) quadrature of the Mellin convolution integral, optionally combined with barycentric interpolation, and a two-dimensional vectorized implementation of the Broda--Kan characteristic-function inversion. The methods are applied to Hake's normalized gain, modeled as a ratio of normal random variables, and are benchmarked against the analytic Hinkley--Marsaglia density evaluated at quadruple precision. The pilot study reports runtimes, speedups, and maximum absolute errors for four representative (a,b) parameter pairs, and the authors provide Jupyter notebooks and code on GitHub.","tokens_in":20718,"tokens_out":16264,"duration_ms":152941,"significance":"The external validation strategy is a clear strength: the reference values come from an independent analytic formula computed at quadruple precision and cross-checked with Arb, rather than from the methods' own outputs, so the reported errors for the 1D DE Mellin-convolution method are meaningful. For the 1D method, the paper convincingly demonstrates that a Numba-based Python DE implementation can match PaCAL's accuracy and approach C-like speeds, and the reproducibility artifacts are a substantial asset. The 2D Broda--Kan implementation addresses a real gap, since no widely available vectorized implementation exists. However, the theoretical justification of the 2D method is not currently secure: the stated inversion theorem requires absolute convergence, a condition that is not verified and is in fact false for the normal cases tested, as discussed below. If that gap is repaired, the paper would be a useful contribution to computational statistics and metrology.","major_comments":[{"comment":"The theorem is invoked as the theoretical basis for the 2D method, but its hypothesis of absolute convergence is not satisfied for the pilot's normal cases. For independent X1 ~ N(a,1) and X2 ~ N(b,1), the integrand in (B4) has a non-removable 1/t singularity at t=0 whose coefficient, Re[phi_{X1}(s)(ib+xs) exp(-ibxs - x^2 s^2 / 2)], is odd in s, so the double integral converges only conditionally after the s-integration cancels the singularity; the integral of the absolute value diverges logarithmically. The manuscript neither verifies absolute convergence nor states an improper/principal-value or iterated-integral interpretation. Consequently, the reported accuracy of the 2D method is not logically connected to the theorem as stated. Please add a rigorous treatment of the conditional convergence (or replace (B4) with an equivalent difference form over t>0 with a removable singularity), and state explicitly which quadrature interpretation Algorithm 1 implements.","section":"Appendix B, Theorem 6 and Eq. (3.13)"},{"comment":"The pseudocode is not self-contained enough to support the reproducibility claim. It does not state the quadrature limits, the step sizes, or the fact that the t2 summation is restricted to positive half-offset points while the t1 summation is symmetric; those choices are exactly what make the conditionally convergent integral in (3.13) well defined. As written, the algorithm is a bare trapezoidal sum of a non-absolutely convergent integral, and a reader cannot determine from the paper alone how N and h were chosen or how convergence was monitored. Please provide the complete quadrature rule and a concrete convergence criterion, and clarify how the singularity cancellation is handled.","section":"Appendix C, Algorithm 1"}],"minor_comments":[{"comment":"The paragraph beginning 'The method's speed, accuracy, and reliability were tested...' appears twice in the conclusions with only minor wording differences; please remove the duplicate.","section":"Section 5"},{"comment":"The sentence 'The only commercial is MATLAB' is grammatically incomplete; it should read 'The only commercial tool is MATLAB.'","section":"Section 4.1"},{"comment":"The software name is spelled both 'PaCAL' and 'PACAL'; please standardize to the package's official capitalization.","section":"Sections 3.1 and 4.2"},{"comment":"The notation in Eq. (2.6) uses delta_2 both for the coefficient of variation of X2 and for the sum delta_1^2 + delta_2^2 in the variance expression; please introduce a distinct symbol for the combined quantity to avoid confusion.","section":"Eq. (2.6)"},{"comment":"In Table 2, the accuracy column reports values such as '3e-16' and '1e-02' without stating that these are maximum absolute errors; please add a note that all accuracy entries are epsilon_max as defined in Section 4.1.","section":"Section 4.1 and Table 2"}],"recommendation":"major_revision","confidential_remarks":"The 1D DE Mellin-convolution contribution is sound and well validated, and the external analytic benchmark is a genuine strength. The main obstacle to publication is the 2D Broda--Kan part: the absolute-convergence hypothesis of Theorem 6 is not met for the tested normal cases, and the algorithm description is too sparse to establish reproducibility on its own. These issues appear fixable within the paper's scope by proving the conditional convergence with an appropriate principal-value/iterated-integral formulation and by supplying a complete algorithm specification. The self-citations to prior DE-quadrature software are tool provenance rather than circular validation."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Two things you should know. First, the 1D double-exponential Mellin convolution with barycentric interpolation is a genuine, well-validated computational result. The pilot benchmark uses an external analytic reference computed at quadruple precision and cross-checked with Arb, and the reported errors support the accuracy claims. Second, the 2D Broda–Kan inversion has a loose thread. Theorem 6 is stated under an absolute-convergence condition, and the normal-normal pilot cases do not satisfy it—the integrand has a 1/t singularity and the double integral is not absolutely convergent. The authors never verify the condition or give a principal-value/conditional-integral interpretation. I don't buy the stronger stress-test claim that Algorithm 1 omits a term and would give half the Cauchy density; for the (a,b)=(0,0) case the single-term B4 integral evaluates to the full Cauchy density at x=0, so the algorithm is probably right. But the theorem as written does not cover the examples used to validate it.\n\nThe novelty is honest. Mellin convolution and Broda–Kan are classical, but the DE-plus-barycentric pairing and the vectorized implementation are not, and the open notebooks are a real service. The 2D method is roughly two orders of magnitude slower than the 1D method, so its practical case rests on generality beyond normal ratios—which the authors themselves label as preliminary. Citations look fair; the self-citations are tool provenance, not a circular load.\n\nThe soft spot that matters is the gap between the stated theorem and the pilot cases. A referee can fix it by adding a short lemma showing the s-integration removes the t=0 singularity, and by being explicit about truncation and conditional convergence. The reproducibility claim would also be stronger with a commit hash or DOI, but that is minor.\n\nWho is this for? Statisticians and metrologists who need fast, accurate ratio PDFs beyond the normal case. I would cite the 1D method; the 2D method is worth watching. The paper deserves a serious referee. I would send it to peer review with a request for revision rather than desk-reject it.","headline":"The 1D Mellin–DE method is the real, well-validated contribution; the 2D Broda–Kan part works empirically but the stated theorem's convergence condition is not satisfied for the pilot cases, so the paper needs a caveat or a repaired argument.","tokens_in":21248,"tokens_out":12334,"would_cite":true,"duration_ms":117056,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["62E15","65D30"],"pacs":[],"model":"deepseek-v4-flash","headline":"The paper claims that two numerical methods—1D double-exponential quadrature of a Mellin convolution and 2D vectorized Broda–Khan characteristic-function inversion—compute ratio-distribution densities with accuracy matching the best…","keywords":["Ratio statistics","Exact probability distributions","Integral transforms","Double exponential quadrature","Hake normalized gain","Effect size","Characteristic function inversion"],"falsifier":"Evaluate the Broda–Khan inversion on a ratio of heavy-tailed variables, for example $X_1\\sim N(0,1)$ and $X_2\\sim t_3$, by comparing the 2D method's output against a high-precision evaluation of the exact Mellin convolution integral. If the error exceeds the reported tolerance or the absolute-convergence condition stated in Theorem 6 is not satisfied, the paper's claim that the 2D method is broadly applicable and reliable would be refuted.","tokens_in":20235,"feed_emoji":"🧮","tokens_out":13359,"duration_ms":117342,"temperature":0.7,"pith_summary":"This paper tries to establish that two numerical techniques can compute the probability density of a ratio of random variables quickly and accurately enough for practical statistics. The first technique applies one-dimensional double-exponential quadrature to the Mellin convolution integral for the ratio's PDF; the second uses a fully vectorized two-dimensional Broda–Khan inversion of characteristic functions. Both are tested on Hake's normalized gain, the standard effect-size measure in physics education, which is the ratio of two correlated normal variables. The authors claim the Python implementation of the first method matches the accuracy of the established PaCAL package and reaches speeds comparable to C code, while the second method supplies the first functional reproducible implementation of Broda–Khan inversion. If these claims hold, researchers can obtain exact ratio-distribution calculations in open tools instead of relying on Monte Carlo or normal approximations.","feed_headline":"Two quadrature methods deliver fast, accurate ratio-distribution PDFs","feed_subtitle":"Both run in open data-science tools and are tested on Hake's normalized gain in effect-size measurement.","key_machinery":"The paper's argument rides on two integral representations of a ratio's density. The first is the Mellin convolution for independent variables, $f_T(t)=\\int_{-\\infty}^{\\infty} f_1(xt)f_2(x)\\,|x|\\,dx$, evaluated by 1D double-exponential quadrature using a $\\sinh$-$\\sinh$ transformation; the DE transformation makes the trapezoidal rule converge exponentially for a broad class of integrands. The second is the Broda–Khan inversion formula, $f_T(x)=\\frac{1}{\\pi^2}\\int_0^\\infty\\int_{-\\infty}^\\infty \\Re\\left[\\varphi_{X_1}(s)\\,\\varphi'_{X_2}(-t-xs)/t\\right]\\,ds\\,dt$, evaluated by a vectorized 2D trapezoidal or Chebyshev quadrature. The Broda–Khan formula is what lets the method work from characteristic functions alone, without requiring known PDFs or independence of the ratio constituents.","core_discovery":"On its own terms, the paper establishes that two numerical routes can replace black-box Monte Carlo for ratio-distribution calculations. The 1D route evaluates the Mellin convolution integral $f_T(t)=\\frac{1}{2\\pi}\\int_{-\\infty}^{\\infty} e^{-((xt-a)^2+(x-b)^2)/2}|x|\\,dx$ with double-exponential quadrature, reaching maximum absolute errors around $10^{-16}$ in sub-millisecond times when compiled with Numba and accelerated with barycentric interpolation. The 2D route evaluates the Broda–Khan inversion integral with a vectorized trapezoidal or Chebyshev quadrature; it is slower by roughly two orders of magnitude but still practically usable, and the authors state it is the first functional reproducible implementation of that inversion. The pilot study on Hake's normalized gain, a ratio of correlated normal variables, shows stable runtimes and errors across four shape regimes, and a preliminary $N(0,1)/\\chi^2(5)$ example indicates the 2D method also works beyond normal constituents. The claim is therefore that these are fast, accurate, and reliable open-tool methods for ratio distributions, not just for this application.","pith_inferences":["If the absolute-convergence condition behind the Broda–Khan formula is checked case by case, the vectorized 2D implementation could give metrology a practical default for uncertainty propagation when the measurement model's characteristic function is known but its density is not.","The 1D DE–Mellin route should transfer directly to other standardized effect-size ratios where the constituents' PDFs are known; a head-to-head run against the same analytic formulas used here would be a cheap validation.","A timing comparison that includes PaCAL's roughly one-second initialization would clarify whether the reported speed parity holds in interactive use, where overhead matters as much as steady-state runtimes.","The pilot study's four parameter pairs cover four qualitative shapes of the ratio distribution, so the natural next step is a systematic parameter sweep to map where each method's quadrature parameters need tuning."],"forward_implications":["For Hake's normalized gain, researchers can now compute exact PDF values numerically in milliseconds instead of relying on Monte Carlo or normal approximations, using the paper's notebooks and code.","The 1D double-exponential Mellin method, when compiled with Numba, reaches speeds comparable to the 20-year-old C code and to PaCAL's Cython implementation, while exposing a single accuracy parameter that trades speed against error.","The 2D Broda–Khan implementation requires only characteristic functions, so it can handle non-normal and negative-valued ratio constituents without assuming independence.","Across the four tested shape regimes of Hake's ratio, both methods kept maximum absolute errors stable and runtimes within a 5% coefficient of variation.","Preliminary results for $X_1\\sim N(0,1)$ over $X_2\\sim\\chi^2(5)$ indicate the Broda–Khan method extends to non-normal ratios, with agreement against a PDF-based approach."],"supporting_citations":[{"why":"supplies the inversion theorem (Appendix B, Thm. 6) that the 2D characteristic-function method implements.","marker":"Broda and Kan (2016)"},{"why":"provides the alternative analytic PDF form and the existing C code that serve as the paper's accuracy and speed baselines.","marker":"Marsaglia (2006)"},{"why":"describes PaCAL, the Mellin-convolution package whose performance the Python DE implementation is benchmarked against.","marker":"(Korze´ n and Jaroszewicz, 2014)"},{"why":"details the bivariate characteristic-function inversion quadrature and vectorization scheme that the 2D implementation adapts.","marker":"Mijanovi´ c et al., 2023"},{"why":"introduces double-exponential quadrature, the numerical engine of the 1D method.","marker":"(Takahasi and Mori, 1974)"},{"why":"provides the analytic PDF and modality results used to cross-check the numerical outputs.","marker":"Pham-Gia et al. (2007)"},{"why":"gives the classical analytic PDF for ratios of correlated normal variables, used as the reference for the Hake-ratio application.","marker":"Hinkley (1969)"},{"why":"offers CharFunTool, whose bivariate inversion routine is the starting point the paper modifies for the 2D Broda–Khan integral.","marker":"Witkovsk´ y (2023)"}],"fun_headline_variants":["Two quadrature methods replace Monte Carlo for ratio PDFs","Fast ratio PDFs via double-exponential and Broda-Khan quadrature","Open-tool quadrature nails accurate ratio distributions","First reproducible Broda-Khan inversion for ratio PDFs"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The second method assumes that the double integral defining the density settles to a finite value before the numerical truncation happens, but the paper never verifies this assumption directly for the cases it runs.","fun_headline_variants_meta":{"raw":{"variants":["Two quadrature methods replace Monte Carlo for ratio PDFs","Fast ratio PDFs via double-exponential and Broda-Khan quadrature","Open-tool quadrature nails accurate ratio distributions","First reproducible Broda-Khan inversion for ratio PDFs"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000335,"raw_usage":{"total_tokens":1887,"prompt_tokens":1005,"completion_tokens":882,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":621,"completion_tokens_details":{"reasoning_tokens":813}},"tokens_in":621,"tokens_out":882,"duration_ms":8419,"temperature":1.0,"reasoning_tokens":813,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-12T17:01:40.488937+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Evaluate the Broda–Khan inversion on a ratio of heavy-tailed variables, for example $X_1\\sim N(0,1)$ and $X_2\\sim t_3$, by comparing the 2D method's output against a high-precision evaluation of the exact Mellin convolution integral. If the error exceeds the reported tolerance or the absolute-convergence condition stated in Theorem 6 is not satisfied, the paper's claim that the 2D method is broadly applicable and reliable would be refuted.","supporting_citations":[],"review_version":1}