From 813b3b173390ea44059a944cdce39a725de28eb2 Mon Sep 17 00:00:00 2001
From: Talen Mudaly <talen.mudaly@student.manchester.ac.uk>
Date: Sun, 12 Jul 2026 22:26:55 +0100
Subject: [PATCH 1/3] Add icdf skeleton for Wald distribution

---
 pymc/distributions/continuous.py | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/pymc/distributions/continuous.py b/pymc/distributions/continuous.py
index d616d72a81..b5877c4ada 100644
--- a/pymc/distributions/continuous.py
+++ b/pymc/distributions/continuous.py
@@ -1072,6 +1072,17 @@ def logcdf(value, mu, lam, alpha):
             msg="mu > 0, lam > 0, alpha >= 0",
         )
 
+    def icdf(value, mu, lam, alpha):
+        res = pt.as_tensor_variable(pt.nan) + value * 0
+        res = check_icdf_value(res, value)
+        return check_icdf_parameters(
+            res,
+            mu > 0,
+            lam > 0,
+            alpha >= 0,
+            msg="mu > 0, lam > 0, alpha >= 0",
+        )
+
 
 class BetaClippedRV(BetaRV):
     @classmethod

From c382f6838691caf66bef07f14ea6440e947d09fd Mon Sep 17 00:00:00 2001
From: Talen Mudaly <talen.mudaly@student.manchester.ac.uk>
Date: Wed, 29 Jul 2026 04:07:26 +0100
Subject: [PATCH 2/3] Add icdf for Wald distribution using Newton's method
 (Giner & Smyth 2016)

---
 pymc/distributions/continuous.py | 30 +++++++++++++++++++++++++++---
 1 file changed, 27 insertions(+), 3 deletions(-)

diff --git a/pymc/distributions/continuous.py b/pymc/distributions/continuous.py
index b5877c4ada..800e84af14 100644
--- a/pymc/distributions/continuous.py
+++ b/pymc/distributions/continuous.py
@@ -1073,10 +1073,34 @@ def logcdf(value, mu, lam, alpha):
         )
 
     def icdf(value, mu, lam, alpha):
-        res = pt.as_tensor_variable(pt.nan) + value * 0
-        res = check_icdf_value(res, value)
+        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(
-            res,
+            result,
             mu > 0,
             lam > 0,
             alpha >= 0,

From c86a83ae5d294531ab12a148081fbd835adefa2d Mon Sep 17 00:00:00 2001
From: Talen Mudaly <talen.mudaly@student.manchester.ac.uk>
Date: Tue, 4 Aug 2026 20:38:28 +0100
Subject: [PATCH 3/3] Add test for Wald icdf

---
 tests/distributions/test_continuous.py | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/tests/distributions/test_continuous.py b/tests/distributions/test_continuous.py
index 9f0c729e2c..85e881bd2c 100644
--- a/tests/distributions/test_continuous.py
+++ b/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",
         [