From 6198c5e08f65a5b93f69f44900ebc25bf80f8999 Mon Sep 17 00:00:00 2001
From: ShreeshaM07 <shreesham2k22@yahoo.com>
Date: Wed, 17 Jul 2024 23:37:22 +0530
Subject: [PATCH 1/8] [ENH] Tweedie Distribution

---
 skpro/distributions/tweedie.py | 152 +++++++++++++++++++++++++++++++++
 1 file changed, 152 insertions(+)
 create mode 100644 skpro/distributions/tweedie.py

diff --git a/skpro/distributions/tweedie.py b/skpro/distributions/tweedie.py
new file mode 100644
index 000000000..baa18aaeb
--- /dev/null
+++ b/skpro/distributions/tweedie.py
@@ -0,0 +1,152 @@
+# copyright: skpro developers, BSD-3-Clause License (see LICENSE file)
+"""Tweedie probability distribution."""
+
+__author__ = ["ShreeshaM07"]
+
+import numpy as np
+
+from skpro.distributions.base import BaseDistribution
+
+# import pandas as pd
+
+
+class Tweedie(BaseDistribution):
+    """Tweedie Distribution.
+
+    Parameters
+    ----------
+    pw : float or array of float (1D or 2D)
+        Power parameter
+    mu : float or array of float (1D or 2D)
+        mean of the normal distribution
+    scale : float or array of float (1D or 2D)
+        scale parameter
+    index : pd.Index, optional, default = RangeIndex
+    columns : pd.Index, optional, default = RangeIndex
+    """
+
+    _tags = {
+        "capabilities:approx": ["pdfnorm"],
+        "capabilities:exact": [
+            "mean",
+            "var",
+            "energy",
+            "pdf",
+            "log_pdf",
+            "cdf",
+            "ppf",
+            "pmf",
+            "log_pmf",
+        ],
+        "distr:measuretype": "mixed",
+        "distr:paramtype": "parametric",
+        "broadcast_init": "on",
+    }
+
+    def __init__(self, pw, mu, scale, index=None, columns=None):
+        # from skpro.distributions.gamma import Gamma
+        from skpro.distributions.normal import Normal
+        from skpro.distributions.poisson import Poisson
+
+        self.pw = pw
+        self.mu = mu
+        self.scale = scale
+
+        if pw == 0:
+            self._norm = Normal(mu=mu, sigma=scale, index=index, columns=columns)
+        elif pw == 1:
+            self._pois = Poisson(mu=mu, index=index, columns=columns)
+        # elif pw == 2:
+        #     self._gam = Gamma(alpha=alpha,beta=beta,index=index,columns=columns)
+
+        super().__init__(index=index, columns=columns)
+
+    def _pdf(self, x):
+        """Probability density function.
+
+        Parameters
+        ----------
+        x : 2D np.ndarray, same shape as ``self``
+            values to evaluate the pdf at
+
+        Returns
+        -------
+        2D np.ndarray, same shape as ``self``
+            pdf values at the given points
+        """
+        from scipy.special import wright_bessel
+
+        pw = self.pw
+        mu = self.mu
+        scale = self.scale
+        if pw == 0:
+            return self._norm.pdf(x)
+        elif pw > 1 and pw < 2:
+            theta = np.power(mu, 1 - pw) / (1 - pw)
+            kappa = np.power(mu, 2 - pw) / (2 - pw)
+            alpha = (2 - pw) / (1 - pw)
+            t = ((pw - 1) * scale / x) ** alpha
+            t /= (2 - pw) * scale
+            a = 1 / x * wright_bessel(-alpha, 0, t)
+            return a * np.exp((x * theta - kappa) / scale)
+
+    def _pmf(self, x):
+        """Probability mass function.
+
+        Private method, to be implemented by subclasses.
+        """
+        pw = self.pw
+        mu = self.mu
+        scale = self.scale
+        if pw == 1:
+            return self._pois.pdf(x)
+        elif pw > 1 and pw < 2:
+            return np.exp(-np.power(mu, 2 - pw) / (scale * (2 - pw)))
+
+    def _log_pdf(self, x):
+        """Logarithmic probability density function.
+
+        Parameters
+        ----------
+        x : 2D np.ndarray, same shape as ``self``
+            values to evaluate the pdf at
+
+        Returns
+        -------
+        2D np.ndarray, same shape as ``self``
+            log pdf values at the given points
+        """
+
+    def _log_pmf(self, x):
+        """Logarithmic probability mass function.
+
+        Private method, to be implemented by subclasses.
+        """
+
+    def _cdf(self, x):
+        """Cumulative distribution function.
+
+        Parameters
+        ----------
+        x : 2D np.ndarray, same shape as ``self``
+            values to evaluate the cdf at
+
+        Returns
+        -------
+        2D np.ndarray, same shape as ``self``
+            cdf values at the given points
+        """
+
+    def _ppf(self, p):
+        """Quantile function = percent point function = inverse cdf.
+
+        Parameters
+        ----------
+        p : 2D np.ndarray, same shape as ``self``
+            values to evaluate the ppf at
+
+        Returns
+        -------
+        2D np.ndarray, same shape as ``self``
+            ppf values at the given points
+        """

From a0eb95db7941770683f2f251de777955367b9dcc Mon Sep 17 00:00:00 2001
From: ShreeshaM07 <shreesham2k22@yahoo.com>
Date: Fri, 19 Jul 2024 01:29:41 +0530
Subject: [PATCH 2/8] compound poisson gamma

---
 skpro/distributions/cmp_pois_gam.py | 130 ++++++++++++++++++++++++++++
 skpro/distributions/tweedie.py      |  98 ++++++++++++++-------
 2 files changed, 199 insertions(+), 29 deletions(-)
 create mode 100644 skpro/distributions/cmp_pois_gam.py

diff --git a/skpro/distributions/cmp_pois_gam.py b/skpro/distributions/cmp_pois_gam.py
new file mode 100644
index 000000000..78d9d8adb
--- /dev/null
+++ b/skpro/distributions/cmp_pois_gam.py
@@ -0,0 +1,130 @@
+# copyright: skpro developers, BSD-3-Clause License (see LICENSE file)
+"""Tweedie probability distribution."""
+
+__author__ = ["ShreeshaM07"]
+
+import numpy as np
+
+from skpro.distributions.base import BaseDistribution
+
+# import pandas as pd
+
+
+class CmpPoissonGamma(BaseDistribution):
+    """Compound Poisson Gamma Distribution.
+
+    Parameters
+    ----------
+    pow : float or array of float (1D or 2D)
+        Power parameter should be in range (1,2)
+    mu : float or array of float (1D or 2D)
+        mean of the normal distribution
+    scale : float or array of float (1D or 2D)
+        scale parameter
+    index : pd.Index, optional, default = RangeIndex
+    columns : pd.Index, optional, default = RangeIndex
+    """
+
+    _tags = {
+        "capabilities:approx": ["pdfnorm"],
+        "capabilities:exact": [
+            "mean",
+            "var",
+            "energy",
+            "pdf",
+            "log_pdf",
+            "cdf",
+            "ppf",
+            "pmf",
+            "log_pmf",
+        ],
+        "distr:measuretype": "mixed",
+        "distr:paramtype": "parametric",
+        "broadcast_init": "on",
+    }
+
+    def __init__(self, pow, mu, scale, index=None, columns=None):
+        self.pow = pow
+        self.mu = mu
+        self.scale = scale
+
+        super().__init__(index=index, columns=columns)
+
+    def _pdf(self, x):
+        """Probability density function.
+
+        Parameters
+        ----------
+        x : 2D np.ndarray, same shape as ``self``
+            values to evaluate the pdf at
+
+        Returns
+        -------
+        2D np.ndarray, same shape as ``self``
+            pdf values at the given points
+        """
+        from scipy.special import wright_bessel
+
+        pow = self.pow
+        mu = self.mu
+        scale = self.scale
+
+        theta = np.power(mu, 1 - pow) / (1 - pow)
+        kappa = np.power(mu, 2 - pow) / (2 - pow)
+        alpha = (2 - pow) / (1 - pow)
+        t = ((pow - 1) * scale / x) ** alpha
+        t /= (2 - pow) * scale
+        a = 1 / x * wright_bessel(-alpha, 0, t)
+        return a * np.exp((x * theta - kappa) / scale)
+
+    def _pmf(self, x):
+        """Probability mass function.
+
+        Private method, to be implemented by subclasses.
+        """
+        pow = self.pow
+        mu = self.mu
+        scale = self.scale
+
+        return np.exp(-np.power(mu, 2 - pow) / (scale * (2 - pow)))
+
+    def _cdf(self, x):
+        """Cumulative distribution function.
+
+        Parameters
+        ----------
+        x : 2D np.ndarray, same shape as ``self``
+            values to evaluate the cdf at
+
+        Returns
+        -------
+        2D np.ndarray, same shape as ``self``
+            cdf values at the given points
+        """
+        from scipy.integrate import quad
+
+        cdf_val = np.array(quad(self._pdf, 0, x)[0])
+        return cdf_val
+
+    def _ppf(self, p):
+        """Quantile function = percent point function = inverse cdf.
+
+        Parameters
+        ----------
+        p : 2D np.ndarray, same shape as ``self``
+            values to evaluate the ppf at
+
+        Returns
+        -------
+        2D np.ndarray, same shape as ``self``
+            ppf values at the given points
+        """
+        from scipy.optimize import brentq
+
+        def objective(x, p):
+            return self._cdf(x) - p
+
+        ppf_val = np.array(
+            brentq(objective, 0, 1000, args=(p,))
+        )  # Adjust the upper bound as necessary
+        return ppf_val
diff --git a/skpro/distributions/tweedie.py b/skpro/distributions/tweedie.py
index baa18aaeb..f42b3c0a8 100644
--- a/skpro/distributions/tweedie.py
+++ b/skpro/distributions/tweedie.py
@@ -15,7 +15,7 @@ class Tweedie(BaseDistribution):
 
     Parameters
     ----------
-    pw : float or array of float (1D or 2D)
+    pow : float or array of float (1D or 2D)
         Power parameter
     mu : float or array of float (1D or 2D)
         mean of the normal distribution
@@ -43,21 +43,29 @@ class Tweedie(BaseDistribution):
         "broadcast_init": "on",
     }
 
-    def __init__(self, pw, mu, scale, index=None, columns=None):
-        # from skpro.distributions.gamma import Gamma
+    def __init__(self, pow, mu, scale, index=None, columns=None):
+        from skpro.distributions.cmp_pois_gam import CmpPoissonGamma
+        from skpro.distributions.gamma import Gamma
         from skpro.distributions.normal import Normal
         from skpro.distributions.poisson import Poisson
 
-        self.pw = pw
+        self.pow = pow
         self.mu = mu
         self.scale = scale
-
-        if pw == 0:
+        mu = np.array(mu)
+        scale = np.array(scale)
+        if pow == 0:
             self._norm = Normal(mu=mu, sigma=scale, index=index, columns=columns)
-        elif pw == 1:
+        elif pow == 1:
             self._pois = Poisson(mu=mu, index=index, columns=columns)
-        # elif pw == 2:
-        #     self._gam = Gamma(alpha=alpha,beta=beta,index=index,columns=columns)
+        elif pow > 1 and pow < 2:
+            self._cmp_pg = CmpPoissonGamma(
+                pow=pow, mu=mu, scale=scale, index=index, columns=columns
+            )
+        elif pow == 2:
+            alpha = (mu / scale) ** 2
+            beta = mu / scale**2
+            self._gam = Gamma(alpha=alpha, beta=beta, index=index, columns=columns)
 
         super().__init__(index=index, columns=columns)
 
@@ -74,34 +82,32 @@ def _pdf(self, x):
         2D np.ndarray, same shape as ``self``
             pdf values at the given points
         """
-        from scipy.special import wright_bessel
+        pow = self.pow
 
-        pw = self.pw
-        mu = self.mu
-        scale = self.scale
-        if pw == 0:
+        if pow == 0:
             return self._norm.pdf(x)
-        elif pw > 1 and pw < 2:
-            theta = np.power(mu, 1 - pw) / (1 - pw)
-            kappa = np.power(mu, 2 - pw) / (2 - pw)
-            alpha = (2 - pw) / (1 - pw)
-            t = ((pw - 1) * scale / x) ** alpha
-            t /= (2 - pw) * scale
-            a = 1 / x * wright_bessel(-alpha, 0, t)
-            return a * np.exp((x * theta - kappa) / scale)
+        elif pow == 1:
+            return self._pois.pdf(x)
+        elif pow > 1 and pow < 2:
+            return self._cmp_pg.pdf(x)
+        elif pow == 2:
+            return self._gam.pdf(x)
 
     def _pmf(self, x):
         """Probability mass function.
 
         Private method, to be implemented by subclasses.
         """
-        pw = self.pw
-        mu = self.mu
-        scale = self.scale
-        if pw == 1:
-            return self._pois.pdf(x)
-        elif pw > 1 and pw < 2:
-            return np.exp(-np.power(mu, 2 - pw) / (scale * (2 - pw)))
+        pow = self.pow
+
+        if pow == 0:
+            return self._norm.pmf(x)
+        elif pow == 1:
+            return self._pois.pmf(x)
+        elif pow > 1 and pow < 2:
+            return self._cmp_pg.pmf(x)
+        elif pow == 2:
+            return self._gam.pmf(x)
 
     def _log_pdf(self, x):
         """Logarithmic probability density function.
@@ -116,12 +122,26 @@ def _log_pdf(self, x):
         2D np.ndarray, same shape as ``self``
             log pdf values at the given points
         """
+        pow = self.pow
+        if pow == 0:
+            return self._norm.log_pdf(x)
+        elif pow == 1:
+            return self._pois.log_pdf(x)
+        elif pow > 1 and pow < 2:
+            return self._cmp_pg.log_pdf(x)
 
     def _log_pmf(self, x):
         """Logarithmic probability mass function.
 
         Private method, to be implemented by subclasses.
         """
+        pow = self.pow
+        if pow == 0:
+            return self._norm.log_pmf(x)
+        elif pow == 1:
+            return self._pois.log_pmf(x)
+        elif pow > 1 and pow < 2:
+            return self._cmp_pg.log_pmf(x)
 
     def _cdf(self, x):
         """Cumulative distribution function.
@@ -136,6 +156,16 @@ def _cdf(self, x):
         2D np.ndarray, same shape as ``self``
             cdf values at the given points
         """
+        pow = self.pow
+
+        if pow == 0:
+            return self._norm.cdf(x)
+        elif pow == 1:
+            return self._pois.cdf(x)
+        elif pow > 1 and pow < 2:
+            return self._cmp_pg.cdf(x)
+        elif pow == 2:
+            return self._gam.cdf(x)
 
     def _ppf(self, p):
         """Quantile function = percent point function = inverse cdf.
@@ -150,3 +180,13 @@ def _ppf(self, p):
         2D np.ndarray, same shape as ``self``
             ppf values at the given points
         """
+        pow = self.pow
+
+        if pow == 0:
+            return self._norm.ppf(p)
+        elif pow == 1:
+            return self._pois.ppf(p)
+        elif pow > 1 and pow < 2:
+            return self._cmp_pg.ppf(p)
+        elif pow == 2:
+            return self._gam.ppf(p)

From 12d9826abc339523d71685e107185dbd4a8c7c3d Mon Sep 17 00:00:00 2001
From: ShreeshaM07 <shreesham2k22@yahoo.com>
Date: Wed, 24 Jul 2024 11:47:31 +0530
Subject: [PATCH 3/8] manual implementation of pdf of cpg

---
 skpro/distributions/cmp_pois_gam.py | 122 ++++++++++++----------------
 1 file changed, 52 insertions(+), 70 deletions(-)

diff --git a/skpro/distributions/cmp_pois_gam.py b/skpro/distributions/cmp_pois_gam.py
index 78d9d8adb..1aad870d8 100644
--- a/skpro/distributions/cmp_pois_gam.py
+++ b/skpro/distributions/cmp_pois_gam.py
@@ -1,26 +1,28 @@
 # copyright: skpro developers, BSD-3-Clause License (see LICENSE file)
-"""Tweedie probability distribution."""
+"""Compound poisson gamma probability distribution."""
 
 __author__ = ["ShreeshaM07"]
 
+import math
+
 import numpy as np
+from scipy.special import gamma as gam_fun
+from scipy.stats import poisson
 
 from skpro.distributions.base import BaseDistribution
 
-# import pandas as pd
-
 
 class CmpPoissonGamma(BaseDistribution):
     """Compound Poisson Gamma Distribution.
 
     Parameters
     ----------
-    pow : float or array of float (1D or 2D)
-        Power parameter should be in range (1,2)
-    mu : float or array of float (1D or 2D)
-        mean of the normal distribution
-    scale : float or array of float (1D or 2D)
-        scale parameter
+    lambda_ : float or array of float (1D or 2D)
+        The rate parameter of the Poisson distribution.
+    alpha : float or array of float (1D or 2D)
+        The shape parameter of the Gamma distribution.
+    beta : float or array of float (1D or 2D)
+        The rate parameter (inverse scale) of the Gamma distribution.
     index : pd.Index, optional, default = RangeIndex
     columns : pd.Index, optional, default = RangeIndex
     """
@@ -43,11 +45,10 @@ class CmpPoissonGamma(BaseDistribution):
         "broadcast_init": "on",
     }
 
-    def __init__(self, pow, mu, scale, index=None, columns=None):
-        self.pow = pow
-        self.mu = mu
-        self.scale = scale
-
+    def __init__(self, lambda_, alpha, beta, index=None, columns=None):
+        self.lambda_ = lambda_
+        self.alpha = alpha
+        self.beta = beta
         super().__init__(index=index, columns=columns)
 
     def _pdf(self, x):
@@ -63,68 +64,49 @@ def _pdf(self, x):
         2D np.ndarray, same shape as ``self``
             pdf values at the given points
         """
-        from scipy.special import wright_bessel
-
-        pow = self.pow
-        mu = self.mu
-        scale = self.scale
-
-        theta = np.power(mu, 1 - pow) / (1 - pow)
-        kappa = np.power(mu, 2 - pow) / (2 - pow)
-        alpha = (2 - pow) / (1 - pow)
-        t = ((pow - 1) * scale / x) ** alpha
-        t /= (2 - pow) * scale
-        a = 1 / x * wright_bessel(-alpha, 0, t)
-        return a * np.exp((x * theta - kappa) / scale)
-
-    def _pmf(self, x):
+        lam = self.lambda_
+        alpha = self.alpha
+        beta = self.beta
+        pdf_value = np.zeros_like(x)
+        tol = 1e-10
+
+        for idx, val in np.ndenumerate(x):
+            if val <= 0:
+                continue  # PDF is zero for non-positive values
+
+            const_term = np.exp(-beta * val) / ((np.exp(lam) - 1) * val)
+            i = 1
+            while True:
+                t1 = lam * (pow(beta * val, alpha))
+                numer = pow(t1, i)
+                i_fact = math.factorial(i)
+                gamma_fun = gam_fun(i * alpha)
+                denom = i_fact * gamma_fun
+
+                term = numer / denom
+                pdf_value[idx] += term
+
+                if term < tol:
+                    break
+                i += 1
+                if i > 1000:  # safeguard to prevent infinite loop
+                    break
+            pdf_value[idx] = pdf_value[idx] * const_term
+
+        return pdf_value
+
+    def _pmf(self, k):
         """Probability mass function.
 
-        Private method, to be implemented by subclasses.
-        """
-        pow = self.pow
-        mu = self.mu
-        scale = self.scale
-
-        return np.exp(-np.power(mu, 2 - pow) / (scale * (2 - pow)))
-
-    def _cdf(self, x):
-        """Cumulative distribution function.
-
         Parameters
         ----------
-        x : 2D np.ndarray, same shape as ``self``
-            values to evaluate the cdf at
+        k : 2D np.ndarray, same shape as ``self``
+            values to evaluate the pmf at
 
         Returns
         -------
         2D np.ndarray, same shape as ``self``
-            cdf values at the given points
+            pmf values at the given points
         """
-        from scipy.integrate import quad
-
-        cdf_val = np.array(quad(self._pdf, 0, x)[0])
-        return cdf_val
-
-    def _ppf(self, p):
-        """Quantile function = percent point function = inverse cdf.
-
-        Parameters
-        ----------
-        p : 2D np.ndarray, same shape as ``self``
-            values to evaluate the ppf at
-
-        Returns
-        -------
-        2D np.ndarray, same shape as ``self``
-            ppf values at the given points
-        """
-        from scipy.optimize import brentq
-
-        def objective(x, p):
-            return self._cdf(x) - p
-
-        ppf_val = np.array(
-            brentq(objective, 0, 1000, args=(p,))
-        )  # Adjust the upper bound as necessary
-        return ppf_val
+        lambda_ = self.lambda_
+        return poisson.pmf(k, lambda_)

From edfa74fe7d3d5e2f991e7a143cee28ab5876707b Mon Sep 17 00:00:00 2001
From: ShreeshaM07 <shreesham2k22@yahoo.com>
Date: Thu, 25 Jul 2024 15:41:31 +0530
Subject: [PATCH 4/8] implement cdf but not giving correct answers

---
 skpro/distributions/cmp_pois_gam.py | 61 ++++++++++++++++++++++++++++-
 1 file changed, 59 insertions(+), 2 deletions(-)

diff --git a/skpro/distributions/cmp_pois_gam.py b/skpro/distributions/cmp_pois_gam.py
index 1aad870d8..62b4c6f57 100644
--- a/skpro/distributions/cmp_pois_gam.py
+++ b/skpro/distributions/cmp_pois_gam.py
@@ -6,8 +6,6 @@
 import math
 
 import numpy as np
-from scipy.special import gamma as gam_fun
-from scipy.stats import poisson
 
 from skpro.distributions.base import BaseDistribution
 
@@ -64,6 +62,8 @@ def _pdf(self, x):
         2D np.ndarray, same shape as ``self``
             pdf values at the given points
         """
+        from scipy.special import gamma as gam_fun
+
         lam = self.lambda_
         alpha = self.alpha
         beta = self.beta
@@ -108,5 +108,62 @@ def _pmf(self, k):
         2D np.ndarray, same shape as ``self``
             pmf values at the given points
         """
+        from scipy.stats import poisson
+
         lambda_ = self.lambda_
         return poisson.pmf(k, lambda_)
+
+    def _compute_crj(self, r, j, rho):
+        from itertools import combinations_with_replacement
+
+        from scipy.special import factorial
+
+        c_rj = 0
+        partitions = [
+            p for p in combinations_with_replacement(range(1, r + 1), j) if sum(p) == r
+        ]
+        for partition in partitions:
+            term = 1
+            for s_i in partition:
+                term *= factorial(rho + 1 + s_i) / (
+                    factorial(rho - 1) * factorial(s_i + 2)
+                )
+            c_rj += term
+        c_rj *= (rho**2 + rho) ** (-r / 2 - j)
+        return c_rj
+
+    def _compute_hr(self, x, r, rho):
+        from scipy.special import eval_hermitenorm, factorial
+
+        hr = np.zeros_like(x, dtype=float)
+        for j in range(1, r + 1):
+            H = eval_hermitenorm(r + 2 * j, x)
+            crj = self._compute_crj(r, j, rho)
+            hr += H * crj / factorial(j)
+        return hr
+
+    def _cdf(self, x):
+        """Cumulative distribution function.
+
+        Parameters
+        ----------
+        x : 2D np.ndarray, same shape as ``self``
+            values to evaluate the cdf at
+
+        Returns
+        -------
+        2D np.ndarray, same shape as ``self``
+            cdf values at the given points
+        """
+        rho = self.alpha
+        max_r = 10  # Adjust as needed
+        phi_x = np.exp(-(x**2) / 2) / np.sqrt(2 * np.pi)
+        series_sum = np.zeros_like(x, dtype=float)
+        for r in range(1, max_r + 1):
+            hr_x = self._compute_hr(x, r, rho)
+            series_sum += hr_x
+            # print(f"r={r}, hr_x={hr_x}, series_sum={series_sum}")
+            # if np.any(np.abs(hr_x) < 1e-6):  # Convergence check
+            #     break
+        cdf = phi_x * (1 + series_sum)
+        return cdf

From daf40d04ec41c8de06f4471a42cb916f7ab44f6b Mon Sep 17 00:00:00 2001
From: ShreeshaM07 <shreesham2k22@yahoo.com>
Date: Tue, 30 Jul 2024 14:00:57 +0530
Subject: [PATCH 5/8] used numpy/scipy functions to speed up cdf

---
 skpro/distributions/cmp_pois_gam.py | 28 +++++++++++++++-------------
 1 file changed, 15 insertions(+), 13 deletions(-)

diff --git a/skpro/distributions/cmp_pois_gam.py b/skpro/distributions/cmp_pois_gam.py
index 62b4c6f57..df79a258e 100644
--- a/skpro/distributions/cmp_pois_gam.py
+++ b/skpro/distributions/cmp_pois_gam.py
@@ -116,20 +116,22 @@ def _pmf(self, k):
     def _compute_crj(self, r, j, rho):
         from itertools import combinations_with_replacement
 
-        from scipy.special import factorial
-
-        c_rj = 0
-        partitions = [
-            p for p in combinations_with_replacement(range(1, r + 1), j) if sum(p) == r
-        ]
-        for partition in partitions:
-            term = 1
-            for s_i in partition:
-                term *= factorial(rho + 1 + s_i) / (
-                    factorial(rho - 1) * factorial(s_i + 2)
-                )
-            c_rj += term
+        from scipy.special import comb
+
+        partitions = np.array(
+            [
+                p
+                for p in combinations_with_replacement(range(1, r + 1), j)
+                if sum(p) == r
+            ]
+        )
+        if partitions.size == 0:
+            return 0
+
+        term = np.prod([comb(rho + 1 + s_i, s_i + 2) for s_i in partitions.T], axis=0)
+        c_rj = np.sum(term)
         c_rj *= (rho**2 + rho) ** (-r / 2 - j)
+
         return c_rj
 
     def _compute_hr(self, x, r, rho):

From d968633ffd8217b018dcd75ea2372f6ea0f50ee9 Mon Sep 17 00:00:00 2001
From: ShreeshaM07 <shreesham2k22@yahoo.com>
Date: Sun, 4 Aug 2024 01:05:08 +0530
Subject: [PATCH 6/8] replaced with normal pdf

---
 skpro/distributions/cmp_pois_gam.py | 10 ++++------
 1 file changed, 4 insertions(+), 6 deletions(-)

diff --git a/skpro/distributions/cmp_pois_gam.py b/skpro/distributions/cmp_pois_gam.py
index df79a258e..15ccfcb73 100644
--- a/skpro/distributions/cmp_pois_gam.py
+++ b/skpro/distributions/cmp_pois_gam.py
@@ -131,7 +131,6 @@ def _compute_crj(self, r, j, rho):
         term = np.prod([comb(rho + 1 + s_i, s_i + 2) for s_i in partitions.T], axis=0)
         c_rj = np.sum(term)
         c_rj *= (rho**2 + rho) ** (-r / 2 - j)
-
         return c_rj
 
     def _compute_hr(self, x, r, rho):
@@ -157,15 +156,14 @@ def _cdf(self, x):
         2D np.ndarray, same shape as ``self``
             cdf values at the given points
         """
+        from scipy.stats import norm
+
         rho = self.alpha
-        max_r = 10  # Adjust as needed
-        phi_x = np.exp(-(x**2) / 2) / np.sqrt(2 * np.pi)
+        max_r = 10
+        phi_x = norm.pdf(x)
         series_sum = np.zeros_like(x, dtype=float)
         for r in range(1, max_r + 1):
             hr_x = self._compute_hr(x, r, rho)
             series_sum += hr_x
-            # print(f"r={r}, hr_x={hr_x}, series_sum={series_sum}")
-            # if np.any(np.abs(hr_x) < 1e-6):  # Convergence check
-            #     break
         cdf = phi_x * (1 + series_sum)
         return cdf

From 9066d6d02ef8242f88cf791672a17dcc84b259c3 Mon Sep 17 00:00:00 2001
From: ShreeshaM07 <shreesham2k22@yahoo.com>
Date: Mon, 12 Aug 2024 23:27:01 +0530
Subject: [PATCH 7/8] Permutation instead of just combination

---
 skpro/distributions/cmp_pois_gam.py | 15 ++++++++++++++-
 1 file changed, 14 insertions(+), 1 deletion(-)

diff --git a/skpro/distributions/cmp_pois_gam.py b/skpro/distributions/cmp_pois_gam.py
index 15ccfcb73..98f3a98af 100644
--- a/skpro/distributions/cmp_pois_gam.py
+++ b/skpro/distributions/cmp_pois_gam.py
@@ -114,7 +114,7 @@ def _pmf(self, k):
         return poisson.pmf(k, lambda_)
 
     def _compute_crj(self, r, j, rho):
-        from itertools import combinations_with_replacement
+        from itertools import combinations_with_replacement, permutations
 
         from scipy.special import comb
 
@@ -125,6 +125,11 @@ def _compute_crj(self, r, j, rho):
                 if sum(p) == r
             ]
         )
+
+        partitions = np.vstack(
+            [np.array(list(set(permutations(sub_list)))) for sub_list in partitions]
+        )
+
         if partitions.size == 0:
             return 0
 
@@ -140,7 +145,15 @@ def _compute_hr(self, x, r, rho):
         for j in range(1, r + 1):
             H = eval_hermitenorm(r + 2 * j, x)
             crj = self._compute_crj(r, j, rho)
+            # print("1. r = ",r," and j =",j)
             hr += H * crj / factorial(j)
+        #     print("2. He(",r+2*j,",",x,") = \n",H)
+        #     print("3. crj for r=",r,"j=",j," = ",crj)
+        #     print("4. j!",factorial(j))
+        #     print("5. h_",r,"upto from j=1 to j=",j," = ",hr)
+        #     print("-------------------------------")
+        # print("6. h_",r," = ",hr)
+        # print()
         return hr
 
     def _cdf(self, x):

From 893dcdff9c96826155fc786bf09c687ecc193997 Mon Sep 17 00:00:00 2001
From: ShreeshaM07 <shreesham2k22@yahoo.com>
Date: Tue, 13 Aug 2024 23:46:29 +0530
Subject: [PATCH 8/8] eqn 3.8 instead of 3.10

---
 skpro/distributions/cmp_pois_gam.py | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/skpro/distributions/cmp_pois_gam.py b/skpro/distributions/cmp_pois_gam.py
index 98f3a98af..6fd44478f 100644
--- a/skpro/distributions/cmp_pois_gam.py
+++ b/skpro/distributions/cmp_pois_gam.py
@@ -143,7 +143,7 @@ def _compute_hr(self, x, r, rho):
 
         hr = np.zeros_like(x, dtype=float)
         for j in range(1, r + 1):
-            H = eval_hermitenorm(r + 2 * j, x)
+            H = eval_hermitenorm(r + 2 * j - 1, x)
             crj = self._compute_crj(r, j, rho)
             # print("1. r = ",r," and j =",j)
             hr += H * crj / factorial(j)
@@ -171,12 +171,14 @@ def _cdf(self, x):
         """
         from scipy.stats import norm
 
+        lambda_ = self.lambda_
         rho = self.alpha
         max_r = 10
-        phi_x = norm.pdf(x)
+        phi_x = np.exp(-(x**2) / 2)
+        Phi_x = norm.cdf(x)
         series_sum = np.zeros_like(x, dtype=float)
         for r in range(1, max_r + 1):
             hr_x = self._compute_hr(x, r, rho)
-            series_sum += hr_x
-        cdf = phi_x * (1 + series_sum)
+            series_sum += pow(lambda_, -r / 2) * hr_x
+        cdf = Phi_x - phi_x * series_sum
         return cdf