Implement icdf for Univariate distribution

pymc/distributions/continuous.py

@@ -1019,8 +1019,41 @@ def logcdf(value, mu, lam, alpha):
             msg="mu > 0, lam > 0, alpha >= 0",
         )
 
-    def icdf(mu, lam):
-        return stats.invgauss(mu=mu, scale=lam)
+    def icdf(value, mu, lam, alpha=0):
+        """Calculate the inverse cumulative distribution function (icdf) of the Wald distribution.
+
+        Args
+        ----
+            value (float): Probability value between 0 and 1.
+            mu (float): Mean of the distribution.
+            lam (float): Scale of the distribution.
+
+        Returns
+        -------
+            float: The value x such that P(W <= x) = p, where W is the Wald distribution with given mu and lam.
+
+        Raises
+        ------
+            ValueError: If value is not between 0 and 1.
+            ValueError: If lam is not positive.
+        """
+        # Compute standard deviation and location parameter
+        std = at.sqrt(lam)
+        loc = alpha + mu
+
+        # Compute inverse standard normal CDF
+        z = at.sqrt(2) * at.erfinv(2 * value - 1)
+
+        # Compute Wald ICDF
+        x = loc + std * z / (1 - 0.5 * std * z)
+
+        return check_parameters(
+            x,
+            0 <= value <= 1,
+            lam > 0,
+            alpha >= 0,
+            msg="0<=val<=1, lam > 0, alpha >= 0",
+        )
 
 
 class BetaClippedRV(BetaRV):