From 59a13dbe7dc294e58fff0e98ee251216d31f809b Mon Sep 17 00:00:00 2001
From: EdoLu <ed.luini@gmail.com>
Date: Sat, 11 Nov 2023 15:31:21 +0100
Subject: [PATCH] distribution moment method compatible with scipy version >=
 1.11.3

---
 gemact/distributions.py | 57 ++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 56 insertions(+), 1 deletion(-)

diff --git a/gemact/distributions.py b/gemact/distributions.py
index 2057848..a213bc1 100644
--- a/gemact/distributions.py
+++ b/gemact/distributions.py
@@ -210,7 +210,11 @@ def moment(self, n):
         """
         hf.assert_type_value(n, 'n', logger, (float, int), lower_bound=1, lower_close=True)
         n = int(n)
-        return self._dist.moment(n=n)
+        try:
+            output = self._dist.moment(n=n)
+        except:
+            output = self._dist.moment(order=n)
+        return output
 
     def skewness(self):
         """
@@ -3341,6 +3345,31 @@ def lev(self, v):
         out = (1 - np.exp(-self.theta * v)) / self.theta
         out[v < 0] = v[v < 0]
         return out
+    
+    def partial_moment(self, n, low, up):
+        """
+        Partial moment of order n.
+
+        :param n: moment order.
+        :type n: ``int``
+        :param low: lower limit of the partial moment.
+        :type low: ``int``, ``float``
+        :param up: upper limit of the partial moment.
+        :type up: ``int``, ``float``
+        :return: raw partial moment of order n.
+        :rtype: ``float`` 
+        """
+
+        hf.assert_type_value(low, 'low', logger, type=(int, float),
+        lower_bound=0, upper_bound=float('inf'), upper_close=False)
+        hf.assert_type_value(up, 'up', logger, type=(int, float), lower_bound=low, lower_close=False)
+        hf.assert_type_value(n, 'n', logger, type=(int, float), lower_bound=0, lower_close=True)
+        n = int(n)
+
+        scale = 1/self.theta
+        output = scale**n * np.exp(special.loggamma(n+1) - special.loggamma(1))
+        output *= stats.gamma(n+1, scale=scale).cdf(up) - stats.gamma(n+1, scale=scale).cdf(low)
+        return output
 
 
 # Gamma
@@ -3431,6 +3460,32 @@ def lev(self, v):
         out[flt] = (alpha / beta) * special.gammainc(alpha + 1, beta * v[flt]) + v[flt] * (1 - special.gammainc(alpha, beta * v[flt]))
         return out
 
+    def partial_moment(self, n, low, up):
+        """
+        Partial moment of order n.
+
+        :param n: moment order.
+        :type n: ``int``
+        :param low: lower limit of the partial moment.
+        :type low: ``int``, ``float``
+        :param up: upper limit of the partial moment.
+        :type up: ``int``, ``float``
+        :return: raw partial moment of order n.
+        :rtype: ``float`` 
+        """
+
+        hf.assert_type_value(low, 'low', logger, type=(int, float),
+        lower_bound=0, upper_bound=float('inf'), upper_close=False)
+        hf.assert_type_value(up, 'up', logger, type=(int, float), lower_bound=low, lower_close=False)
+        hf.assert_type_value(n, 'n', logger, type=(int, float), lower_bound=0, lower_close=True)
+        n = int(n)
+
+        scale = self.scale
+        shape = self.a
+        output = scale**n * np.exp(special.loggamma(shape + n) - special.loggamma(shape))
+        output *= stats.gamma(shape + n, scale=scale).cdf(up) - stats.gamma(shape + n, scale=scale).cdf(low)
+        return output
+
 
 # Inverse Gamma
 class InvGamma(_ContinuousDistribution):