Add icdf for Wald distribution

pymc/distributions/continuous.py

@@ -1072,6 +1072,41 @@ def logcdf(value, mu, lam, alpha):
             msg="mu > 0, lam > 0, alpha >= 0",
         )
 
+    def icdf(value, mu, lam, alpha):
+        phi = lam / mu
+        kappa = 3 * mu / (2 * lam)
+        m = mu * (pt.sqrt(1 + kappa**2) - kappa)
+        res = pt.switch(
+            pt.and_(pt.gt(value, 1e-5), pt.lt(value, 1 - 1e-5)),
+            m,
+            pt.switch(
+                pt.gt(value, 1 - 1e-5),
+                Gamma.icdf(value=value, alpha=phi, scale=mu / phi),
+                mu / (phi * Normal.icdf(value=value, mu=0, sigma=1) ** 2),
+            ),
+        )
+
+        def newton_step(q, *_):
+            cdf_q = pt.exp(Wald.logcdf(q, mu, lam, 0))
+            pdf_q = pt.exp(Wald.logp(q, mu, lam, 0))
+            return q + (value - cdf_q) / pdf_q
+
+        q, _ = pytensor.scan(
+            newton_step,
+            outputs_info=[res],
+            n_steps=20,
+        )
+
+        result = q[-1] + alpha
+        result = check_icdf_value(result, value)
+        return check_icdf_parameters(
+            result,
+            mu > 0,
+            lam > 0,
+            alpha >= 0,
+            msg="mu > 0, lam > 0, alpha >= 0",
+        )
+
 
 class BetaClippedRV(BetaRV):
     @classmethod

tests/distributions/test_continuous.py

@@ -348,6 +348,14 @@ def test_wald_logcdf(self):
             lambda value, mu, alpha: st.invgauss.logcdf(value, mu=mu, loc=alpha),
         )
 
+    def test_wald_icdf(self):
+        check_icdf(
+            pm.Wald,
+            {"mu": Rplus, "alpha": Rplus},  # no lam — defaults to 1
+            lambda q, mu, alpha: st.invgauss.ppf(q, mu=mu, loc=alpha),
+            decimal=4,
+        )
+
     @pytest.mark.parametrize(
         "value,mu,lam,phi,alpha,logp",
         [