From 8822346958cdf7586dad054dc874468231e4d784 Mon Sep 17 00:00:00 2001
From: Param Thakkar <128291516+ParamThakkar123@users.noreply.github.com>
Date: Wed, 10 Jan 2024 09:14:22 +0530
Subject: [PATCH 1/3] Added icdf for Rice and skewnormal distributions

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

diff --git a/pymc/distributions/continuous.py b/pymc/distributions/continuous.py
index 76932db6f6..3374d4e15d 100644
--- a/pymc/distributions/continuous.py
+++ b/pymc/distributions/continuous.py
@@ -76,6 +76,9 @@ def polyagamma_cdf(*args, **kwargs):
 from scipy import stats
 from scipy.interpolate import InterpolatedUnivariateSpline
 from scipy.special import expit
+from scipy.stats import norm
+from scipy.optimize import newton
+from scipy.special import ive, iv
 
 from pymc.distributions import transforms
 from pymc.distributions.dist_math import (
@@ -2993,6 +2996,23 @@ def logp(value, mu, sigma, alpha):
             tau > 0,
             msg="tau > 0",
         )
+    def icdf(prob, mu, sigma, alpha, max_iter=100, tol=1e-8):
+        def cdf_difference(x):
+            return norm.cdf(x) - 2 * norm.cdf(alpha * x) - prob
+        
+        x0 = norm.ppf(prob)
+        x = x0
+        
+        for _ in range(max_iter):
+            diff = cdf_difference(x)
+            
+            if np.abs(diff) < tol:
+                return mu + sigma * x
+            
+            derivative = norm.pdf(x) - 2 * alpha * norm.pdf(alpha * x)
+            x -= derivative
+        
+        return mu + sigma * x
 
 
 class Triangular(BoundedContinuous):
@@ -3348,6 +3368,17 @@ def logp(value, b, sigma):
             b >= 0,
             msg="sigma >= 0, b >= 0",
         )
+    def icdf(prob, nu, sigma, x0=1.0, max_iter=100, tol=1e-8):
+        def cdf(x):
+            return (1 + x / sigma**2) * np.exp(-x**2 / (2 * sigma**2)) * iv(0, x * nu / sigma**2) - prob
+        
+        def cdf_derivative(x):
+            return ((1 - x**2 / sigma**2) * np.exp(-x**2 / (2 * sigma**2)) * iv(0, x * nu / sigma**2)
+                + (x / sigma**2) * np.exp(-x**2 / (2 * sigma**2)) * ive(0, x * nu / sigma**2)) * nu / sigma**2
+        
+        approx_icdf = newton(cdf, x0, fprime=cdf_derivative, tol=tol, maxiter=max_iter)
+
+        return approx_icdf
 
 
 class Logistic(Continuous):

From 8baf68bca996430274bfe535ff3960dc22630404 Mon Sep 17 00:00:00 2001
From: Param Thakkar <128291516+ParamThakkar123@users.noreply.github.com>
Date: Fri, 12 Jan 2024 18:53:52 +0530
Subject: [PATCH 2/3] Update continuous.py

---
 pymc/distributions/continuous.py | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/pymc/distributions/continuous.py b/pymc/distributions/continuous.py
index 3374d4e15d..a8c0f148a6 100644
--- a/pymc/distributions/continuous.py
+++ b/pymc/distributions/continuous.py
@@ -3006,7 +3006,7 @@ def cdf_difference(x):
         for _ in range(max_iter):
             diff = cdf_difference(x)
             
-            if np.abs(diff) < tol:
+            if pt.abs(diff) < tol:
                 return mu + sigma * x
             
             derivative = norm.pdf(x) - 2 * alpha * norm.pdf(alpha * x)
@@ -3370,17 +3370,16 @@ def logp(value, b, sigma):
         )
     def icdf(prob, nu, sigma, x0=1.0, max_iter=100, tol=1e-8):
         def cdf(x):
-            return (1 + x / sigma**2) * np.exp(-x**2 / (2 * sigma**2)) * iv(0, x * nu / sigma**2) - prob
+            return (1 + x / sigma**2) * pt.exp(-x**2 / (2 * sigma**2)) * iv(0, x * nu / sigma**2) - prob
         
         def cdf_derivative(x):
-            return ((1 - x**2 / sigma**2) * np.exp(-x**2 / (2 * sigma**2)) * iv(0, x * nu / sigma**2)
-                + (x / sigma**2) * np.exp(-x**2 / (2 * sigma**2)) * ive(0, x * nu / sigma**2)) * nu / sigma**2
+            return ((1 - x**2 / sigma**2) * pt.exp(-x**2 / (2 * sigma**2)) * iv(0, x * nu / sigma**2)
+                + (x / sigma**2) * pt.exp(-x**2 / (2 * sigma**2)) * ive(0, x * nu / sigma**2)) * nu / sigma**2
         
         approx_icdf = newton(cdf, x0, fprime=cdf_derivative, tol=tol, maxiter=max_iter)
 
         return approx_icdf
 
-
 class Logistic(Continuous):
     r"""
     Logistic log-likelihood.

From 5a2265ffd7287e60c163b91c55c16c4c7cf8efba Mon Sep 17 00:00:00 2001
From: Param Thakkar <128291516+ParamThakkar123@users.noreply.github.com>
Date: Fri, 12 Jan 2024 19:05:20 +0530
Subject: [PATCH 3/3] Update continuous.py

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

diff --git a/pymc/distributions/continuous.py b/pymc/distributions/continuous.py
index a8c0f148a6..46aa8f0836 100644
--- a/pymc/distributions/continuous.py
+++ b/pymc/distributions/continuous.py
@@ -3012,7 +3012,13 @@ def cdf_difference(x):
             derivative = norm.pdf(x) - 2 * alpha * norm.pdf(alpha * x)
             x -= derivative
         
-        return mu + sigma * x
+        res = mu + sigma * x
+        res = check_icdf_value(res, prob)
+        return check_parameters(
+            res,
+            sigma > 0,
+            msg="sigma > 0"
+        )
 
 
 class Triangular(BoundedContinuous):
@@ -3377,8 +3383,13 @@ def cdf_derivative(x):
                 + (x / sigma**2) * pt.exp(-x**2 / (2 * sigma**2)) * ive(0, x * nu / sigma**2)) * nu / sigma**2
         
         approx_icdf = newton(cdf, x0, fprime=cdf_derivative, tol=tol, maxiter=max_iter)
-
-        return approx_icdf
+        res = approx_icdf
+        res = check_icdf_value(res, prob)
+        return check_parameters(
+            res,
+            sigma > 0,
+            msg="sigma > 0"
+        )
 
 class Logistic(Continuous):
     r"""