FEAT: add skew t distribution

gwpopulation/utils.py

@@ -183,6 +183,42 @@ def logsubexp(log_p, log_q):
     return xp.nan_to_num(xp.exp(log_pdf)) * (xx >= low) * (xx <= high)
 
 
+
+def skewt(xx, aa, bb, loc=0, scale=1):
+    r"""
+    Jones and Faddy skew-t distribution (implementation based on :code:`scipy`).
+
+    .. math::
+
+        z &= \frac{x - \mu}{\sigma} \\
+        p(x) &= \frac{1}{\sigma C_{a,b}} \left(1 + \frac{z}{\sqrt{a + b + z^2}}\right)^{a + \frac{1}{2}} \left(1 - \frac{z}{\sqrt{a + b + z^2}}\right)^{b + \frac{1}{2}} \\
+        C_{a,b} &= {2^{a + b - 1} B(a, b) \sqrt{a + b}}
+
+    Parameters
+    ----------
+    xx: float, array-like
+        The abscissa values (:math:`x`)
+    aa: float
+        The first shape parameter (:math:`a`)
+    bb: float
+        The second shape parameter (:math:`b`)
+    loc: float
+        The location parameter (:math:`\mu`)
+    scale: float
+        The scale parameter (:math:`\sigma`)
+
+    Returns
+    -------
+    prob: float, array-like
+        The distribution evaluated at `xx`
+    """
+    zz = (xx - loc) / scale
+    c = 2 ** (aa + bb - 1) * scs.beta(aa, bb) * (aa + bb)**0.5
+    d1 = (1 + zz / (aa + bb + zz ** 2)**0.5) ** (aa + 0.5)
+    d2 = (1 - zz / (aa + bb + zz ** 2)**0.5) ** (bb + 0.5)
+    return d1 * d2 / c / scale
+
+
 def unnormalized_2d_gaussian(xx, yy, mu_x, mu_y, sigma_x, sigma_y, covariance):
     r"""
     Compute the probability distribution for a correlated 2-dimensional Gaussian

test/utils_test.py

@@ -127,7 +127,27 @@ def test_truncnorm_matches_scipy(backend):
         assert max(abs(gwpop_vals - scipy_vals)) < 1e-3
 
 
-def test_matches_scipy(backend):
+def test_skewt_matches_scipy(backend):
+    from scipy.stats import fj_skew_t
+
+    gwpopulation.set_backend(backend)
+    xp = gwpopulation.utils.xp
+    xx = xp.linspace(-2, 2, 1000)
+    for ii in range(N_TEST):
+        mu = np.random.uniform(-10, 10)
+        sigma = np.random.uniform(0, 5)
+        aa = np.random.uniform(0, 100)
+        bb = np.random.uniform(0, 100)
+        gwpop_vals = utils.to_numpy(
+            utils.skewt(xx, aa=aa, bb=bb, loc=mu, scale=sigma)
+        )
+        scipy_vals = fj_skew_t(
+            loc=mu, scale=sigma, a=aa, b=bb
+        ).pdf(utils.to_numpy(xx))
+        assert max(abs(gwpop_vals - scipy_vals)) < 1e-3
+
+
+def test_vonmises_matches_scipy(backend):
     gwpopulation.set_backend(backend)
     xp = gwpopulation.utils.xp
     xx = xp.linspace(0, 2 * np.pi, 1000)