Fix: Incorrect use of partial in TweedieDistribution._rowwise_gradient_hessian

CHANGELOG.rst

@@ -10,6 +10,11 @@ Changelog
 UNRELEASED
 ----------
 
+**Bug fix:
+
+- Fixed a bug where :meth:`~glum.TweedieDistribution._rowwise_gradient_hessian` and :meth:`~glum.TweedieDistribution._eta_mu_deviance` would call functions with wrong arguments in the p = 3 case.
+- Fixed :class:`glum.InverseGaussianDistribution` not using the optimized gradient, hessian and deviance implementations, as well as those gradients and hessians being incorrect.
+
 **Other changes:**
 
 - Build and test with Python 3.13 in CI.

src/glum/_distribution.py

@@ -671,8 +671,8 @@ def _rowwise_gradient_hessian(
             f = gamma_log_rowwise_gradient_hessian
         elif 1 < self.power < 2 and isinstance(link, LogLink):
             f = partial(tweedie_log_rowwise_gradient_hessian, p=self.power)
-        elif self.power == 3:
-            f = partial(inv_gaussian_log_rowwise_gradient_hessian, p=self.power)
+        elif self.power == 3 and isinstance(link, LogLink):
+            f = inv_gaussian_log_rowwise_gradient_hessian
 
         if f is not None:
             return f(y, sample_weight, eta, mu, gradient_rows, hessian_rows)
@@ -703,7 +703,7 @@ def _eta_mu_deviance(
         elif 1 < self.power < 2 and isinstance(link, LogLink):
             f = partial(tweedie_log_eta_mu_deviance, p=self.power)
         elif self.power == 3 and isinstance(link, LogLink):
-            f = partial(inv_gaussian_log_eta_mu_deviance, p=self.power)
+            f = inv_gaussian_log_eta_mu_deviance
 
         if f is not None:
             return f(cur_eta, X_dot_d, y, sample_weight, eta_out, mu_out, factor)
@@ -1153,6 +1153,11 @@ def unit_deviance(self, y, mu):  # noqa D
     def _rowwise_gradient_hessian(
         self, link, y, sample_weight, eta, mu, gradient_rows, hessian_rows
     ):
+        if isinstance(link, LogLink):
+            return inv_gaussian_log_rowwise_gradient_hessian(
+                y, sample_weight, eta, mu, gradient_rows, hessian_rows
+            )
+
         return super()._rowwise_gradient_hessian(
             link, y, sample_weight, eta, mu, gradient_rows, hessian_rows
         )
@@ -1169,8 +1174,8 @@ def _eta_mu_deviance(
         mu_out,
     ):
         if isinstance(link, LogLink):
-            return tweedie_log_eta_mu_deviance(
-                cur_eta, X_dot_d, y, sample_weight, eta_out, mu_out, factor, p=3.0
+            return inv_gaussian_log_eta_mu_deviance(
+                cur_eta, X_dot_d, y, sample_weight, eta_out, mu_out, factor
             )
 
         return super()._eta_mu_deviance(

src/glum/_functions.pyx

@@ -276,8 +276,8 @@ def inv_gaussian_log_rowwise_gradient_hessian(
         inv_mu = 1 / mu[i]
         inv_mu2 = inv_mu ** 2
 
-        gradient_rows_out[i] = 2 * weights[i] * (inv_mu - y[i] * inv_mu2)
-        hessian_rows_out[i] = 2 * weights[i] * (inv_mu - 2 * y[i] * inv_mu2)
+        gradient_rows_out[i] = weights[i] * (y[i] * inv_mu2 - inv_mu)
+        hessian_rows_out[i] = weights[i] * inv_mu
 
 def inv_gaussian_log_likelihood(
     const_floating1d y,

tests/glm/test_distribution.py

@@ -170,6 +170,7 @@ def test_deviance_zero(family, chk_values):
         (InverseGaussianDistribution(), LogLink()),
         (TweedieDistribution(power=1.5), LogLink()),
         (TweedieDistribution(power=2.5), LogLink()),
+        (TweedieDistribution(power=3), LogLink()),
         (BinomialDistribution(), LogitLink()),
         (TweedieDistribution(power=1.5), TweedieLink(1.5)),
         (TweedieDistribution(power=2.5), TweedieLink(2.5)),
@@ -228,6 +229,7 @@ def f(coef2):
         (GammaDistribution(), LogLink(), True),
         (InverseGaussianDistribution(), LogLink(), False),
         (TweedieDistribution(power=1.5), LogLink(), True),
+        (TweedieDistribution(power=3), LogLink(), False),
         (TweedieDistribution(power=4.5), LogLink(), False),
         (NegativeBinomialDistribution(theta=1.0), LogLink(), False),
     ],