REVIEW 2 major objections 5 minor 36 references
Novel computational approaches for ratio distributions with an application to Hake's ratio in effect size measurement
T0 review · 2 major / 5 minor · reviewed 2026-08-12 · deepseek-v4-flash
Pith's one-line read 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…
desk verdict 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. 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 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.
What would settle it
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.
Extended reading notes
Core claim
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.
Load-bearing premise
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.
Editorial extensions
If this is right
- 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.
Reading between the lines
- 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.
Signed reviews
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
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.
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 (2)
- [Appendix B, Theorem 6 and Eq. (3.13)] 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.
- [Appendix C, Algorithm 1] 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.
minor comments (5)
- [Section 5] 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 4.1] The sentence 'The only commercial is MATLAB' is grammatically incomplete; it should read 'The only commercial tool is MATLAB.'
- [Sections 3.1 and 4.2] The software name is spelled both 'PaCAL' and 'PACAL'; please standardize to the package's official capitalization.
- [Eq. (2.6)] 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 4.1 and Table 2] 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.
Circularity Check
No significant circularity: numerical methods are validated against an independent analytic formula; self-citations are tool provenance, not load-bearing.
full rationale
The paper's central claims are that its 1D DE Mellin-convolution implementation and its 2D vectorized Broda-Kan inversion compute ratio PDFs accurately and quickly. Neither claim reduces to its own inputs. Accuracy is assessed by comparing against the external closed-form Hinkley/Marsaglia/Pham-Gia formula (2.3) with (2.5), computed independently in SageMath with quadruple precision and cross-checked against Arb; no parameter of the proposed methods is fitted to the target PDF. The pilot cases (a,b) are fixed inputs, not fitted outputs. Self-citations (Gajdos et al. 2021, Hancova et al. 2022) are used only to attribute the provenance of DE-quadrature code and prior usage, not as evidence for the numerical benchmark results. The Broda-Kan inversion theorem is cited to an external source (Broda and Kan 2016). The possible convergence-condition concerns about Theorem 6 and Appendix C are correctness risks about the validity of the integral representation, not circularity: they do not make the derivation equivalent to its inputs. Accordingly no circular step is identified, and the score is 0.
Assumptions & free parameters
assumptions (6)
- standard math Mellin convolution formula for the PDF of a ratio of independent random variables (Theorem 3, Eq. 3.8).
- standard math Closed-form PDF of the ratio of bivariate normal random variables in terms of Kummer 1F1 or erf (Theorems 1 and 2, Eqs. 2.2 to 2.5).
- standard math Gil-Pelaez inversion formulas for PDF and CDF from a characteristic function (Eq. 3.10).
- domain assumption Broda-Kan ratio inversion theorem requires finite mean, absolutely integrable CF, no atom at X - rY = 0, and absolute convergence of the integrals (Theorem 6, Appendix B, Eqs. B3 and B4).
- domain assumption Scores in Hake's ratio are normal or approximately normal via the central limit theorem (Proposition 4 and Remark 8, Appendix B).
- standard math Double-exponential quadrature converges exponentially for the integrands considered, with optimality results as in Sugihara (1997).
Cite this review
Pith. "Pith review of Novel computational approaches for ratio distributions with an application to Hake's ratio in effect size measurement." pith.science (2026). https://pith.science/paper/OUIQV77Z
@misc{pith2026241112938,
author = {Pith},
title = {Pith review of: Novel computational approaches for ratio distributions with an application to Hake's ratio in effect size measurement},
year = {2026},
howpublished = {\url{https://pith.science/paper/OUIQV77Z}},
note = {Machine review of arXiv:2411.12938}
}
read the original abstract
Ratio statistics and distributions are fundamental in various disciplines, including linear regression, metrology, nuclear physics, operations research, econometrics, biostatistics, genetics, and engineering. In this work, we introduce two novel computational approaches for evaluating ratio distributions using open data science tools and modern numerical quadratures. The first approach employs 1D double exponential quadrature of the Mellin convolution with/without barycentric interpolation, which is a very fast and efficient quadrature technique. The second approach utilizes 2D vectorized Broda-Khan numerical inversion of characteristic functions. It offers broader applicability by not requiring knowledge of PDFs or the independence of ratio constituents. The pilot numerical study, conducted in the context of Hake's ratio - a widely used measure of effect size and educational effectiveness in physics education - demonstrates the proposed methods' speed, accuracy, and reliability. The analytical and numerical explorations also provide more clarifying insight into the theoretical and empirical properties of Hake's ratio distribution. The proposed methods appear promising in a robust framework for fast and exact ratio distribution computations beyond normal random variables, with potential applications in multidimensional statistics and uncertainty analysis in metrology, where precise and reliable data handling is essential.
Reference graph
Works this paper leans on
-
[6]
Phys Rev Phys Educ Res 16 (1): 010108
Why normalized gain should continue to be used in analyzing preinstruction and postinstruction scores on concept inventories. Phys Rev Phys Educ Res 16 (1): 010108. https://doi.org/10.1103/PhysRevPhysEducRes. 16.010108. 22 Cumming, G. and R. Calin-Jageman
-
[10]
Testing independence under a block compound symmetry covariance structure. Stat Papers 64 (2): 677–704. https: //doi.org/10.1007/s00362-022-01335-7. Gajdoˇ s, A., J. Hanˇ c, and M. Hanˇ cov´ a
-
[15]
J Stat Comput Simul 92 (11): 2205–2232
A practical, effective calculation of gamma difference distributions with open data science tools. J Stat Comput Simul 92 (11): 2205–2232. https://doi.org/10.1080/00949655.2021.2023873. 23 Harris, C.R. et al
arXiv 2021
-
[16]
https://doi.org/10.1038/s41586-020-2649-2
Array programming with NumPy.Nature 585 (7825): 357–362. https://doi.org/10.1038/s41586-020-2649-2. Hinkley, D.V
-
[20]
ACM Commun Comput Algebra 47 (3/4): 166–169
Arb: a C library for ball arithmetic. ACM Commun Comput Algebra 47 (3/4): 166–169. https://doi.org/10.1145/2576802.2576828. Johansson, F
-
[24]
https: //doi.org/10.1119/10.0002062
The Physics Teacher 58 (7): 465–469. https: //doi.org/10.1119/10.0002062. Mijanovi´ c, A., B.V. Popovi´ c, and V. Witkovsk´ y
-
[25]
A numerical inversion of the bivariate characteristic function. Appl Math Comput 443: 127807. https://doi.org/ 24 10.1016/j.amc.2022.127807. Mori, M
arXiv 2022
-
[27]
J Franklin Inst 348 (4): 810–822
On the linear combination, product and ratio of normal and Laplace random variables. J Franklin Inst 348 (4): 810–822. https: //doi.org/10.1016/j.jfranklin.2011.01.005. Nissen, J.M. and et al
Show all 36 references
-
[32]
Appl Math Comput 240: 122–139
Kullback–Leibler life time testing. Appl Math Comput 240: 122–139. https://doi.org/10.1016/j.amc.2014. 04.027. Stein, W.A. and others
2014 doi
-
[37]
Stat Probab Lett 96: 223–231
Logarithmic Lambert W×F random variables for the family of chi-squared distributions and their applications. Stat Probab Lett 96: 223–231. https://doi.org/10.1016/j.spl.2014.09.028. 26
2014 doi
- [346]
- [683]
-
[1932]
Biometrika 24 (3-4): 428–440
The distribution of the index in a normal bivariate population. Biometrika 24 (3-4): 428–440. https://doi.org/10.1093/biomet/24.3-4.428. Filipiak, K., M. John, and D. Klein
-
[1951]
Biometrika 38 (3-4): 481–482
Note on the inversion theorem. Biometrika 38 (3-4): 481–482. https://doi.org/10.1093/biomet/38.3-4.481. Gradshteyn, I.S. and I.M. Ryzhik
-
[1965]
J Am Stat Assoc 60 (309): 193–204
Ratios of Normal Variables and Ratios of Sums of Uniform Vari- ables. J Am Stat Assoc 60 (309): 193–204. https://doi.org/10.1080/01621459.1965. 10480783. Marsaglia, G
1965
-
[1969]
Biometrika 56 (3): 635–639
On the Ratio of Two Correlated Normal Random Variables. Biometrika 56 (3): 635–639. https://doi.org/10.2307/2334671. Hinkley, D.V
-
[1973]
Biometrika 60 (2): 415–417
Numerical inversion of a characteristic function. Biometrika 60 (2): 415–417. https://doi.org/10.1093/biomet/60.2.415. D ´ ıaz-Franc´ es, E. and F.J. Rubio
-
[1974]
Publ Res Inst Math Sci 9 (3): 721–741
Double exponential formulas for numerical inte- gration. Publ Res Inst Math Sci 9 (3): 721–741. https://doi.org/10.2977/prims/ 1195192451. Van Rossum, G. and F.L. Drake
-
[1997]
Numer Math 75 (3): 379–395
Optimality of the double exponential formula - functional analysis approach. Numer Math 75 (3): 379–395. https://doi.org/10.1007/s002110050244. 25 Takahasi, H. and M. Mori
-
[1998]
American Journal of Physics 66 (1): 64–74
Interactive-engagement versus traditional methods: A six-thousand- student survey of mechanics test data for introductory physics courses. American Journal of Physics 66 (1): 64–74. https://doi.org/10.1119/1.18809. Hake, R.R
-
[2006]
J Stat Softw 16 (4): 1–10
Ratios of Normal Variables. J Stat Softw 16 (4): 1–10. https: //doi.org/10.18637/jss.v016.i04. MathWorks
-
[2007]
Commun Stat Theory Methods 35 (9): 1569–1591
Density of the Ratio of Two Normal Random Variables and Applications. Commun Stat Theory Methods 35 (9): 1569–1591. https://doi.org/10.1080/03610920600683689. Renyi, A
-
[2010]
London: SAGE Publications
Encyclopedia of Research Design. London: SAGE Publications. Savinainen, A. and P. Scott. 2002, January. Using the Force Concept Inventory to monitor student learning and to plan teaching. Phys. Educ. 37 (1): 53–58. https: //doi.org/10.1088/0031-9120/37/1/307. Severini, T.A
2002 doi
-
[2011]
Comput Sci Eng 13 (2): 31–39
Cython: The best of both worlds. Comput Sci Eng 13 (2): 31–39. https://doi.org/10.1109/MCSE.2010.118. Bowman, N.D
2010 doi
-
[2012]
SIAM J Sci Comput 34 (3): A1241–A1265
Arithmetic Operations on Independent Random Variables: A Numerical Approach. SIAM J Sci Comput 34 (3): A1241–A1265. https: //doi.org/10.1137/110839680. Johansson, F
-
[2013]
Stat Papers 54 (2): 309–323
On the existence of a normal approximation to the distribution of the ratio of two independent normal random variables. Stat Papers 54 (2): 309–323. https://doi.org/10.1007/s00362-012-0429-2. Eaton, J.W
-
[2014]
Qual Reliab Eng Int 30 (8): 1361–1377
Statistical Performance of a Control Chart for Individual Observations Monitoring the Ratio of Two Normal Variables. Qual Reliab Eng Int 30 (8): 1361–1377. https://doi.org/10.1002/qre.1558. Cohen, J
-
[2015]
https://doi
What might psychologists learn from Scholarship of Teaching and Learning in physics? Scholarsh Teach Learn Psychol 1 (1): 100–106. https://doi. org/10.1037/stl0000022. Hanˇ c, J., D. Borovsk´ y, and M. Hanˇ cov´ a
-
[2016]
Acta IMEKO 5 (3): 32–44
Numerical inversion of a characteristic function: An alternative tool to form the probability distribution of output quantity in linear measurement models. Acta IMEKO 5 (3): 32–44. https://doi.org/10.21014/ACTA IMEKO.V5I3
-
[2017]
Communication Research Reports 34 (3): 187–190
The Importance of Effect Size Reporting in Communication Research Reports. Communication Research Reports 34 (3): 187–190. https://doi. org/10.1080/08824096.2017.1353338. Broda, S.A. and R. Kan
2017
-
[2018]
Phys Rev Phys Educ Res 14 (1): 010115
Comparison of normalized gain and Cohen’s d for analyzing gains on concept inventories. Phys Rev Phys Educ Res 14 (1): 010115. https://doi.org/10.1103/PhysRevPhysEducRes.14.010115. Ooura, T
-
[2019]
ACM Trans Math Softw 45 (3): 30:1–30:26
Computing Hypergeometric Functions Rigorously. ACM Trans Math Softw 45 (3): 30:1–30:26. https://doi.org/10.1145/3328732. Korze´ n, M. and S. Jaroszewicz
-
[2020]
Phys Rev Phys Educ Res 16 (1): 010141
Examination of quantitative meth- ods for analyzing data from concept inventories. Phys Rev Phys Educ Res 16 (1): 010141. https://doi.org/10.1103/PhysRevPhysEducRes.16.010141. Celano, G., P. Castagliola, A. Faraz, and S. Fichera
-
[2021]
Environ Res Lett 17 (2): 024011
Environmental concentrations as ratios of random variables. Environ Res Lett 17 (2): 024011. https://doi.org/10.1088/1748-9326/ ac4a9f. Pham-Gia, T., N. Turkkan, and E. Marchand
-
[2022]
https://doi.org/10.1109/ICETA57911.2022.9974868
Interactive Jupyter Notebooks with SageMath in Number Theory, Algebra, Calculus, and Numerical Methods.2022 Pro- ceedings of ICETA: 178–183. https://doi.org/10.1109/ICETA57911.2022.9974868. Galassi, M., J. Davies, J. Theiler, B. Gough, and G. Jungman
2022
-
[2023]
Phys Rev Phys Educ Res 19 (1): 010111
Evidence for a normal distribution of normalized gains. Phys Rev Phys Educ Res 19 (1): 010111. https://doi.org/10.1103/PhysRevPhysEducRes.19. 010111. Coletta, V.P. and J.J. Steinert
Reviewed August 12, 2026 · model on record in the stance chip above.
Discussion (0). Continue with ORCID to comment.