Added random for frank

copulae/archimedean/frank.py

@@ -5,6 +5,8 @@
 from copulae.special.debye import debye_1, debye_2
 from copulae.special.optimize import find_root
 from copulae.special.special_func import log1mexp, log1pexp, poly_log
+from copulae.stats import random_uniform
+from copulae.stats.log import random_log_series_ln1p
 from copulae.types import Array
 from copulae.utility import array_io, as_array
 from .abstract import AbstractArchimedeanCopula
@@ -141,6 +143,7 @@ def pdf(self, u: Array, log=False):
     def psi(self, s):
         assert not np.isnan(self.params), "Copula must have parameters to calculate psi"
 
+        s = np.asarray(s)
         if self.params <= -36:
             return -log1pexp(-s - self.params) / self.params
         elif self.params < 0:
@@ -149,12 +152,30 @@ def psi(self, s):
             return np.exp(-s)
         else:
             const = log1mexp(self.params)
-            if s < const:
-                return np.nan
-            return -log1mexp(s - log1mexp(self.params)) / self.params
+            m = np.less(s, const, where=~np.isnan(s))
+
+            s[m] = np.nan
+            s[~m] = -log1mexp(s[~m] - log1mexp(self.params)) / self.params
+            return s.item(0) if s.size == 1 else s
 
     def random(self, n: int, seed: int = None):
-        pass
+        u = random_uniform(n, self.dim, seed)
+        if abs(self.params) < 1e-7:
+            return u
+
+        if self.dim == 2:
+            v = u[:, 1]
+            a = -abs(self.params)
+            v = -1 / a * np.log1p(-v * np.expm1(-a) / (np.exp(-a * u[:, 0]) * (v - 1) - v))
+            u[:, 1] = 1 - v if self.params > 0 else v
+            return u
+
+        # alpha too large
+        if log1mexp(self.params) == 0:
+            return np.ones((n, self.dim))
+
+        fr = random_log_series_ln1p(-self.params, n)[:, None]
+        return self.psi(-np.log(u) / fr)
 
     @property
     def rho(self):

copulae/stats/log.py

@@ -0,0 +1,45 @@
+import numpy as np
+
+from copulae.types import Numeric
+
+
+def random_log_series_ln1p(alpha: float, size: Numeric = 1):
+    """"""
+    if alpha < -3:
+        # random variate following the logarithmic (series) distribution Log(alpha)
+        p = - alpha / np.log(1 - np.expm1(alpha))
+        u = np.random.uniform(size=size)
+
+        up = u > p
+        x = 1
+        while np.any(up):
+            u[up] -= p
+            x += 1
+            p *= alpha * (x - 1) / x
+
+            up = u > p
+        return u
+    else:
+        # random variate following the logarithmic (series) distribution Log(1 - e^h)
+        e = -np.expm1(alpha)
+        u1 = np.random.uniform(size=size)
+        u2 = np.random.uniform(size=size)
+
+        h = u2 * alpha
+        q = -np.expm1(h)
+        log_q = np.where(h > -0.69314718055994530942, np.log(-np.expm1(h)), np.log1p(-np.exp(h)))
+
+        mask1 = u1 > e
+        mask2 = u1 < q * q
+        mask3 = mask2 & (log_q == 0)
+        mask4 = mask2 & (log_q != 0)
+        mask5 = ~mask2 & (u1 > q)
+        mask6 = ~mask2 & (u1 <= q)
+
+        u1[mask6] = 2
+        with np.errstate(invalid='ignore'):
+            u1[mask4] = 1 + np.floor(np.log(u1 / log_q))[mask4]
+        u1[mask3] = np.inf
+        u1[mask1 | mask5] = 1
+
+        return u1

tests/archimedean/test_clayton.py

@@ -159,6 +159,10 @@ def test_psi(copula, U5):
     assert_array_almost_equal(psi, target, 3)
 
 
+def test_gaussian_random_generates_correctly(copula):
+    assert copula.random(5000).shape == (5000, copula.dim)
+
+
 def test_summary(fitted_clayton):
     smry = fitted_clayton.summary()
     assert isinstance(str(smry), str)

tests/archimedean/test_frank.py

@@ -120,6 +120,10 @@ def test_psi(theta, expected):
         assert np.isnan(copula.psi(log1mexp(theta) - 1e-6))
 
 
+def test_gaussian_random_generates_correctly(copula):
+    assert copula.random(10).shape == (10, copula.dim)
+
+
 @pytest.mark.parametrize('theta, expected', [
     (1.295, 0.211156879535562),
     (1e-9, 1.66666666666667e-10)

tests/archimedean/test_gumbel.py

@@ -122,6 +122,10 @@ def test_gumbel_pdf(copula, U5):
     assert_array_almost_equal(log_cdf, np.log(expected_cdf), DP)
 
 
+def test_gaussian_random_generates_correctly(copula):
+    assert copula.random(10).shape == (10, copula.dim)
+
+
 def test_summary(fitted_gumbel):
     smry = fitted_gumbel.summary()
     assert isinstance(str(smry), str)