From 9aa042d830359e7248ddc5affd34ff3d283430ea Mon Sep 17 00:00:00 2001
From: Malte Londschien <61679398+mlondschien@users.noreply.github.com>
Date: Fri, 7 Jun 2024 16:55:20 +0200
Subject: [PATCH] Use numba_scipy to speed up approximation of CLR critical
 value functio (#80)

* Use numba_scipy to speed up approximation of CLR critical value function.

* numba-scipy not numba_scipy.

* Without numba

* Forgot some lineprof decorators.
---
 ivmodels/tests/conditional_likelihood_ratio.py | 17 +++++++++++------
 1 file changed, 11 insertions(+), 6 deletions(-)

diff --git a/ivmodels/tests/conditional_likelihood_ratio.py b/ivmodels/tests/conditional_likelihood_ratio.py
index c0f7e33..a4a3a26 100644
--- a/ivmodels/tests/conditional_likelihood_ratio.py
+++ b/ivmodels/tests/conditional_likelihood_ratio.py
@@ -82,20 +82,25 @@ def conditional_likelihood_ratio_critical_value_function(
         return 1 - scipy.stats.chi2(q).cdf(z)
 
     if method in ["numerical_integration"]:
-        beta = scipy.stats.beta((q - p) / 2, p / 2)
-        chi2 = scipy.stats.chi2(q)
+        alpha = (q - p) / 2.0
+        beta = p / 2.0
+        k = q / 2.0
+        z_over_2 = z / 2.0
+
         a = s_min / (z + s_min)
+        const = np.power(a, -alpha - beta + 1) / scipy.special.beta(alpha, beta)
 
-        def integrand(b):
-            return beta.pdf(b) * chi2.cdf(z / (1 - a * b))
+        def integrand(y):
+            return const * scipy.special.gammainc(k, z_over_2 / y)
 
         res = scipy.integrate.quad(
             integrand,
-            0,
+            1 - a,
             1,
+            weight="alg",
+            wvar=(beta - 1, alpha - 1),
             epsabs=tol,
         )
-
         return 1 - res[0]
 
     elif method == "power_series":