Probability distributions and calculations for Hake’s ratio statistics in measuring effect size

When Hake introduced the normalized gain in his landmark meta-analysis article [24] to measure instruction quality, he also used empirical data from 62 introductory physics courses ($N=6542$, average $n\approx 100$ participants per course), as shown in the plot on the left in Fig. 1. This visualization also provided an important geometric interpretation with some consequences about used teaching methods.

At point $P$, representing the results of a particular course, the Hake ratio is the tangent of the angle formed by the line connecting point $P=(\overline{S}_{\text{pre}},\overline{S}_{\text{post}}-\overline{S}_{\text{pre}})$ with point $M=(100,0)$. It was observed that courses taught with traditional teaching methods achieved results up to $\hat{g}=0.3$, concentrated around a slope of $\hat{g}_{trad}=0.23$, while courses with interactive teaching methods fell within the range between lines with slopes from 0.3 to 0.6, centered around a slope of $\hat{g}_{int}=0.48$.

Hake’s statistic was and is understood almost exclusively as an empirical measure, representing the practical relative improvement of a group of students after instruction, with recurring discussions about the appropriateness and usefulness of its application [25, 9].

From a statistical perspective, however, only a few articles address the statistical properties of Hake’s statistic and its theoretical or empirical validation. For example, in [8], the author provides empirical evidence showing that real $\hat{g}$ data does not reject the hypothesis of a normal distribution. Regarding theoretical statistical considerations, to the best of our knowledge, only one article [5] attempts such an approach, but a significant weakness of this derivation is its lack of alignment with existing statistical literature.

Therefore, our main goal is to examine the statistical properties and probabilistic computation methods for the Hake normalized gain, offering comprehensive, mutually supportive analytical, empirical, computational, and numerical perspectives to deepen understanding and insights into this ratio in real data applications.

First of all, as stated in the following proposition, using elementary statistical properties of distributions and their moments, we can easily prove that Hake’s ratio $\hat{g}$ is a ratio of two correlated normal random variables (RVs).

In [28], we can find the analytic forms for the PDF and CDF of the ratio $W$ (the CDF of the ratio $W$, expressed via the bivariate normal integral, is again in Appendix A (Theorem 6).

Combining the results of Theorems 1 and 3 provides us the following corollary with the formulas for $a$ and $b$ related to Hake’s ratio.

Another surprising and statistically inconvenient property of the distribution of the ratio of two normal random variables $X_{1}$ and $X_{2}$ is its heavy-tailed nature and lack of finite moments. However, in the case of a unimodal shape, practical applications like Hake’s ratio [8] result in a distribution very close to a normal distribution.

When the means are positive, using a Taylor series expansion of the ratio around the point $(\mu_{1},\mu_{2})$, the first-order term provides the mean approximation, and the second-order term approximates the variance [12, 39]. These "pseudo-moments" are given in the following expression (corresponding to the law of uncertainty propagation):

The conditions for such a reasonable approximation can be formulated numerically [39] and by the existence theorem [12]. This work also introduces a graphical tool, the ROC curve, which is useful for assessing the closeness of the distribution of a specific ratio to the proposed normal approximation.

The specific normal distribution approximation of Hake’s empirical distribution [8] can be seen in Figure 3.

As long-time, experienced users of R and Python for nearly every statistical purpose, we naturally expected an open-source package for calculating the PDF, CDF, approximations, or performing statistical inference related to the ratio of $W$ or $T$. The third surprise for us was discovering that such a specialized package does not exist. The only available software is a nearly 20-year-old $C$ code from the previously mentioned study [39].

The choice of computational method and an appropriate digital tool for scientific computing is therefore left to the potential user. However, when dealing with the analytic formula — containing the hypergeometric confluent function — various computational issues may arise [33, 26]. Therefore, we decided to investigate more deeply the potential numerical approaches that would be suitable for ratio distributions.

As the first result of our research, we found an open-source solution capable of handling probabilistic and statistical computations for our ratio distribution. Specifically, a more general Python package, PaCAL (the Probabilistic CALculator) [35], is able to numerically perform arithmetic operations ($+,-,*,/$) on independent continuous random variables, running on Cython (a C-optimized version of Python). In exploring its documentation, we found that PaCAL computes the ratio distribution through the Mellin transform theory, using the Clenshaw-Curtis (CC) quadrature with barycentric interpolation (via Chebyshev polynomials) [31].

The Mellin transform, defined as [52]:

Arguably the most comprehensive handbook on the topic of the Mellin transform is Springer’s book [52]. The most comprehensive tables of properties and closed expressions for several integral transforms, including the Mellin transform, can be found in the Caltech Project, led by editor A. Erdélyi [13, 14].

It is worth noting that the Mellin transform exists for a real function $f(x)$ that is defined and single-valued almost everywhere for $x\geqslant 0$ and is absolutely integrable over the range $(0,\infty)$. The Mellin transform can be extended [52] to the range $(-\infty,\infty)$. The importance of the Mellin transform in studying the product and ratio of independent random variables lies in the fact that they can be expressed using the so-called Mellin convolution.

In practice, however, the integral can be computed numerically using an appropriate numerical method, quadrature. The mentioned Cython PACAL package calculates such an integral using the CC quadrature with barycentric interpolation.

The probability distribution of any real-valued RV is completely determined by its characteristic function $\varphi_{X}(t)=E\left\{e^{itX}\right\}$, which is mathematically equivalent to the Fourier transform [51]. In many situations, the form of the characteristic function (CF) is analytically known or much simpler than the PDF or CDF of the RV.

In such cases, the characteristic function provides an alternative numerical method for effectively computing the PDF and CDF of the RV through the conceptually simple approach known as the numerical inversion of the characteristic function, which remains less widely known but has been periodically highlighted in research over the past 50 years (e.g., [11, 58, 60]). This computational approach is based on the numerical quadrature of the Gil-Pelaez inversion formulae [21]:

The result of our research and theoretical considerations in Sections 3.1 and 3.2 is that we now have two fundamentally suitable approaches to numerically compute the ratio $X_{1}/X_{2}$ of two independent continuous RVs $X_{1}$ and $X_{2}$. Both approaches are applicable under very mild conditions, i.e., without the necessity of assuming that the RVs to be normally distributed. The first approach is the 1D Mellin convolution integral (3.9) using the PDFs of $X_{1}$ and $X_{2}$. The second approach relies on knowing the CFs of $X_{1}$ and $X_{2}$, involving the calculation of the 2D numerical CF inversion integral (3.14).

From a computational perspective, both integral transforms can be evaluated using any suitable numerical quadrature method. The open source package PaCAL [35] employs the very effective CC quadrature with barycentric interpolation based on Chebyshev polynomials for the numerical calculation of (3.9) [31].

We found that 2D numerical CF inversion integrals can be computed numerically using free MTB’s CharFunTool by [61] with the bivariate algorithm cf2Dist2D [see theoretical details in 42]. This algorithm also employs CC quadrature with barycentric interpolation in 2D and additionally offers a simple 2D trapezoidal rule (TR) quadrature as an alternative option.

It is important to note that the output of cf2Dist2D is the PDF of the bivariate distribution of the random vector $(X_{1},X_{2})$, calculated from its bivariate CF $\phi$ by

To compute our 2D integral (3.14), it is sufficiently to choose $f_{T}(x)=f(x,1)$, with the bivariate characteristic function $\phi(s,t)=2e^{i(x_{1}t_{1}+x_{2}t_{2})}\varphi_{X_{1}}(t_{1})\varphi_{X_{2}}^{\prime}(-x_{1}t_{1}-x_{2}t_{2})/t_{2}.$ The derivative of $\varphi_{X_{2}}$ can be calculated symbolically or numerically. Key details of adjusting the authors’ numerical algorithm for calculating the integral (3.14) are provided in the Appendix (Numerical Quadrature Algorithm).

Regarding numerical quadrature, a very promising alternative to the CC algorithm is the double exponential (DE) quadrature [56, 43]. We remind only the core idea (for more details, see Remark 4), which lies in applying the simple trapezoidal rule to an integral transformed by a suitable DE transformation ($x=\Phi(t);\ a,b\in\mathcal{R}\cup\left\{-\infty,+\infty\right\}$):

so the trapezoidal rule is modified by introducing a weight function $w=\Phi^{\prime}(t)$ with a double-exponential decay rate $\left|\Phi^{\prime}(t)\right|\underset{t\rightarrow\pm\infty}{\rightarrow}\exp\left(-\beta\cdot\exp(\gamma|t|)\right).$ Such DE transformations (three examples are in Table 1) make the trapezoidal rule optimally efficient, yielding exponentially fast and accurate results for a very broad class of functions (even with multiple singularities).

In the Scientific Python ecosystem, the DE quadrature can be seamlessly combined with barycentric interpolation-based Chebyshev polynomials (via the BarycentricInterpolator function from scipy). For open-source computational implementations, the sinh-sinh DE transformation is particularly suitable for the 1D Mellin convolution integral (3.9) in our case. This approach utilizes our Python implementation of the 1D DE quadrature [18], which is a modified Python version of the original C implementation by [47].

Additionally, we have developed a small Python package, Chebyshev.py, specifically tailored for our PDF calculations, enabling the combination of DE quadrature with barycentric interpolation. Regarding the 2D numerical CF inversion (3.14), we decided to adapt the Witkovský’s algorithmic approach for numerical inversion of bivariate CF. As a potential alternative, we also decided to explore MPFUN2020 library [1, 59], a thread-safe, arbitrary-precision Fortran package. This library includes a 2D implementation of the DE quadrature and supports its parallel execution on CPU/GPUs.

For purposes of numerical study benchmarking chosen computational methods and digital tools, we use the open Jupyter interactive computing environment, combined with open-source tools like the computer algebra system SageMath as its kernel, enabling to work with multiple programming languages, such as Python and R, within a unified digital workspace. Additionally, systems like SageMath provide the computational power of traditional software without financial barriers, fostering accessibility and empowering research without the constraints of commercial software. A detailed list of the open digital tools used in our study is in Appendix B.

All computations were conducted on a PC with Windows 11 and a 12th Gen Intel(R) Core(TM) i7-12700K processor @ 3.60 GHz, featuring Intel(R) UHD Graphics 770, 8 performance cores, 20 threads, and 32 GB of RAM. The Python, Fortran, SageMath, and other free open software were installed from their official repositories, while MATLAB (MTB) computations were performed using the version R2023b as a kernel in our Jupyter environment [40].

We conducted one numerical experiment for each available computational method and corresponding digital tool. Each experiment involved 3 runs, with each run containing 10 realizations of computing $T$ at 1000 points uniformly distributed over a chosen interval. This setup resulted in a total of $3\times 10^{4}$ calculations per experiment. Built-in commands or functions were used to measure the execution time of each experiment. To verify accuracy, we generated results using SageMath with quadruple precision (34 digits of accuracy) and cross-checked them by the arbitrary-precision C library Arb.

The benchmark results from each experiment include average runtime, runtime standard deviation, acceleration (relative to the SageMath analytic approach), and accuracy (maximum absolute error). The complete numerical study encompasses systematically carried-out 251 experiments (250 successful and 1 unsuccessful), with some chosen results summarized in Table 2, Figures 5, and 5. All resulting data and code, including 20 Jupyter notebooks, describing also our developments of algorithms in Python and MTB, are available at https://github.com/JupyterPER/HakeRatio.

The benchmark for specific methods and tools consists of:

• •

  13 experiments using the analytic expression (2.6) in SageMath via the fast_float routine (53-bit, $\varepsilon_{\text{rel}}=10^{-3}$ to $10^{-15}$);

• •

  26 experiments using CC quadrature of the 1D Mellin convolution integral (3.9) in PaCAL, with and without Cython initialization (53-bit, $\varepsilon_{\text{rel}}=10^{-3}$ to $10^{-15}$);

• •

  78 experiments using DE quadrature of the 1D Mellin convolution integral (3.9) in Python, with and without Numba, parallelization, barycentric interpolation (53-bit, $\varepsilon_{\text{rel}}=10^{-3}$ to $10^{-15}$);

• •

  16 experiments using trapezoidal rule (TR) of the 2D numerical CF inversion integral (3.14) in Python (53-bit) dealing with different algorithms for an approximating integral sum (compiled in different Numba settings – with/without vectorization, parallization, symbolic/numerical derivative);

• •

  118 experiments using TR or CC quadrature (without/without vectorization, parallelization, symbolic/numerical derivative and various numbers of Chebyshev nodes; using built-in integrators and/or analytical implementations) of the 2D numerical CF inversion integral (3.14) in MTB (53-bit);

• •

  1 experiment using DE quadrature of the 2D numerical CF inversion integral (3.14) in Fortran via MPFUN (53-bit).

In Table 2, we present a sample of deliberately selected experiments divided into three main groups based on the method of computing the $T$ ratio, each group clearly demonstrating the key results of our numerical study. Some of these results are also visualized in Figures 5 and 5.

As for the analytic formula (2.6), the built-in implementations of Kummer’s confluent hypergeometric function for computing the PDF of the $T$ ratio are highly accurate and reliable, which is not always the case for hypergeometric functions [33, 26]. Notably, the SageMath version is an order of magnitude faster (about 7x) than MTB, despite MTB’s reputation for numerical efficiency.

For the Mellin convolution integral, we conducted successful experiments in Python using the PACAL package and Python-based DE implementations. Both CC and DE quadrature methods, combined with high-performance compilers (Cython, Numba), achieved comparable speeds. While PACAL was the fastest, this excludes the 1-second initialization of the barycentric interpolant. Both quadratures delivered reliable results, with DE offering efficient accuracy control, where increasing accuracy by 12 orders of magnitude slows computation by less than one order of magnitude.

At first glance, our newly implemented DE quadrature combined with barycentric interpolation of adjustable accuracy appears to be several times slower. However, it is important to note that the built-in scipy implementation of the barycentric interpolator does not utilize process parallelization. Without parallelization, the speed of the DE quadrature alone is an order of magnitude lower. Ultimately, the use of Chebyshev interpolation accelerated the computation.

The majority of our numerical experiments were conducted using custom MTB and Python algorithms for the numerical inversion of the ratio via the Broda-Khan approach (3.14). While this 2D method is more computationally intensive, achieving speeds 2 orders of magnitude lower than the 1D Mellin convolution, the accuracy and performance are still practically acceptable. At more reduced accuracy, performance is comparable to 1D Python’s DE quadrature or MTB’s analytical computations.

In Python, parallelizing three nested loops (TR, Numba-par) improved 2D inversion speeds by 5-7x per loop, proportional to the number of cores. Without parallelization, the speed drops by 3 orders of magnitude. In MTB, full vectorization (TR, par vect) reduced computation time significantly (from 7.97 s to 0.19 s), proving much more efficient than MTB’s built-in parallelization. With barycentric interpolation, MTB is fast comparably to its analytical ${}_{1}F_{1}$ computations.

One set of experiments—the 2D-DE quadrature in the Fortran library MPFUN2020—was unsuccessful, with inconsistent results and calculation failures.¹¹1We attribute this not to the failure of method but to our limited familiarity with Fortran. Future solutions include collaborating with experts, improving our understanding, or rewriting the code in more familiar languages, potentially with AI assistance, as we did for Marsaglia’s C code.