Repository: skypyproject/skypy
Pull request #210 (2020-07-07 14:39:26)

simplified Smail redshift sampling
## Description

Implement the Smail redshift distribution as a simple function that samples redshifts. This removes the dependency on scipy which throws warnings with v1.5+.

## Checklist
- [x] Follow the [Contributor Guidelines](https://github.com/skypyproject/skypy/blob/master/CONTRIBUTING.md)
- [x] Write unit tests
- [x] Write documentation strings
- [x] Assign someone from your working team to review this pull request
- [x] Assign someone from the infrastructure team to review this pull request


rrjbca: @ntessore can you update your branch and then this is ready to merge

Patch:
diff --git a/skypy/galaxy/redshift.py b/skypy/galaxy/redshift.py
index ba3116f..8da23c6 100644
--- a/skypy/galaxy/redshift.py
+++ b/skypy/galaxy/redshift.py
@@ -5,7 +5,6 @@ models.
 """
 
 import numpy as np
-from scipy import stats
 import scipy.special as sc
 import scipy.integrate
 
@@ -19,7 +18,7 @@ __all__ = [
 ]
 
 
-class smail_gen(stats.rv_continuous):
+def smail(z_median, alpha, beta, size=None):
     r'''Redshifts following the Smail et al. (1994) model.
 
     The redshift follows the Smail et al. [1]_ redshift distribution.
@@ -32,6 +31,9 @@ class smail_gen(stats.rv_continuous):
         Power law exponent (z/z0)^\alpha, must be positive.
     beta : float or array_like of floats
         Log-power law exponent exp[-(z/z0)^\beta], must be positive.
+    size : None or int or tuple
+        Size of the output. If `None`, the size is inferred from the arguments.
+        Default is None.
 
     Notes
     -----
@@ -43,7 +45,7 @@ class smail_gen(stats.rv_continuous):
         p(z) \sim \left(\frac{z}{z_0}\right)^\alpha
                     \exp\left[-\left(\frac{z}{z_0}\right)^\beta\right] \;.
 
-    This is the generalised gamma distribution `scipy.stats.gengamma`.
+    This is the generalised gamma distribution.
 
     References
     ----------
@@ -52,64 +54,18 @@ class smail_gen(stats.rv_continuous):
 
     Examples
     --------
-    >>> from skypy.galaxy.redshift import smail
-
     Sample 10 random variates from the Smail model with `alpha = 1.5` and
     `beta = 2` and median redshift `z_median = 1.2`.
 
-    >>> redshift = smail.rvs(1.2, 1.5, 2.0, size=10)
-
-    Fix distribution parameters for repeated use.
+    >>> from skypy.galaxy.redshift import smail
+    >>> redshift = smail(1.2, 1.5, 2.0, size=10)
 
-    >>> redshift_dist = smail(1.2, 1.5, 2.0)
-    >>> redshift_dist.median()
-    1.2
-    >>> redshift = redshift_dist.rvs(size=10)
     '''
 
-    def _rvs(self, zm, a, b):
-        sz, rs = self._size, self._random_state
-        k = (a+1)/b
-        t = zm**b/sc.gammainccinv(k, 0.5)
-        g = stats.gamma.rvs(k, scale=t, size=sz, random_state=rs)
-        return g**(1/b)
-
-    def _pdf(self, z, zm, a, b):
-        return np.exp(self._logpdf(z, zm, a, b))
-
-    def _logpdf(self, z, zm, a, b):
-        k = (a+1)/b
-        z0 = zm/sc.gammainccinv(k, 0.5)**(1/b)
-        lognorm = np.log(b) - np.log(z0) - sc.gammaln(k)
-        return lognorm + sc.xlogy(a, z/z0) - (z/z0)**b
-
-    def _cdf(self, z, zm, a, b):
-        k = (a+1)/b
-        t = sc.gammainccinv(k, 0.5)
-        return sc.gammainc(k, t*(z/zm)**b)
-
-    def _ppf(self, q, zm, a, b):
-        k = (a+1)/b
-        t = sc.gammainccinv(k, 0.5)
-        return zm*(sc.gammaincinv(k, q)/t)**(1/b)
-
-    def _sf(self, z, zm, a, b):
-        k = (a+1)/b
-        t = sc.gammainccinv(k, 0.5)
-        return sc.gammaincc(k, t*(z/zm)**b)
-
-    def _isf(self, q, zm, a, b):
-        k = (a+1)/b
-        t = sc.gammainccinv(k, 0.5)
-        return zm*(sc.gammainccinv(k, q)/t)**(1/b)
-
-    def _munp(self, n, zm, a, b):
-        k = (a+1)/b
-        z0 = zm/sc.gammainccinv(k, 0.5)**(1/b)
-        return z0**n*np.exp(sc.gammaln((a+n+1)/b) - sc.gammaln(k))
-
-
-smail = smail_gen(a=0., name='smail', shapes='z_median, alpha, beta')
+    k = (alpha+1)/beta
+    t = z_median**beta/sc.gammainccinv(k, 0.5)
+    g = np.random.gamma(shape=k, scale=t, size=size)
+    return g**(1/beta)
 
 
 def herbel_redshift(alpha, a_phi, b_phi, a_m, b_m, cosmology, low=0.0,


Test patch:
diff --git a/skypy/galaxy/tests/test_redshift.py b/skypy/galaxy/tests/test_redshift.py
index 0320cf2..894650d 100644
--- a/skypy/galaxy/tests/test_redshift.py
+++ b/skypy/galaxy/tests/test_redshift.py
@@ -1,76 +1,35 @@
 import numpy as np
-from scipy import stats
-from scipy.stats.tests.common_tests import (
-    check_normalization, check_moment, check_mean_expect, check_var_expect,
-    check_skew_expect, check_kurt_expect, check_edge_support,
-    check_random_state_property, check_pickling)
 
 from skypy.galaxy.redshift import smail, herbel_redshift, herbel_pdf
 from unittest.mock import patch
-import scipy.stats
+from scipy.stats import kstest, norm
 import scipy.integrate
 
 
 def test_smail():
-    # freeze a distribution with parameters
-    args = (1.3, 2.0, 1.5)
-    dist = smail(*args)
-
-    # check that PDF is normalised
-    check_normalization(smail, args, 'smail')
-
-    # check CDF and SF
-    assert np.isclose(dist.cdf(3.) + dist.sf(3.), 1.)
-
-    # check inverse CDF and SF
-    assert np.isclose(dist.ppf(dist.cdf(4.)), 4.)
-    assert np.isclose(dist.isf(dist.sf(5.)), 5.)
-
-    # check median matches parameter
-    zm = np.random.rand(10)
-    assert np.allclose(smail.median(zm, 2.0, 1.5), zm)
-
-    # check moments
-    m, v, s, k = dist.stats(moments='mvsk')
-    check_mean_expect(smail, args, m, 'smail')
-    check_var_expect(smail, args, m, v, 'smail')
-    check_skew_expect(smail, args, m, v, s, 'smail')
-    check_kurt_expect(smail, args, m, v, k, 'smail')
-    check_moment(smail, args, m, v, 'smail')
-
-    # check other properties
-    check_edge_support(smail, args)
-    check_random_state_property(smail, args)
-    check_pickling(smail, args)
-
     # sample a single redshift
-    rvs = dist.rvs()
+    rvs = smail(1.3, 2.0, 1.5)
     assert np.isscalar(rvs)
 
     # sample 10 reshifts
-    rvs = dist.rvs(size=10)
+    rvs = smail(1.3, 2.0, 1.5, size=10)
     assert rvs.shape == (10,)
 
     # sample with implicit sizes
     zm, a, b = np.ones(5), np.ones((7, 5)), np.ones((13, 7, 5))
-    rvs = smail.rvs(z_median=zm, alpha=a, beta=b)
+    rvs = smail(z_median=zm, alpha=a, beta=b)
     assert rvs.shape == np.broadcast(zm, a, b).shape
 
-    # check sampling against own CDF
-    D, p = stats.kstest(smail.rvs, smail.cdf, args=args, N=1000)
-    assert p > 0.01, 'D = {}, p = {}'.format(D, p)
-
     # check sampling, for alpha=0, beta=1, the distribution is exponential
-    D, p = stats.kstest(smail.rvs(0.69315, 1e-100, 1., size=1000), 'expon')
+    D, p = kstest(smail(0.69315, 1e-100, 1., size=1000), 'expon')
     assert p > 0.01, 'D = {}, p = {}'.format(D, p)
 
     # check sampling, for beta=1, the distribution matches a gamma distribution
-    D, p = stats.kstest(smail.rvs(2.674, 2, 1, size=1000), 'gamma', args=(3,))
+    D, p = kstest(smail(2.674, 2, 1, size=1000), 'gamma', args=(3,))
     assert p > 0.01, 'D = {}, p = {}'.format(D, p)
 
     # check sampling, the distribution is a generalised gamma distribution
-    D, p = stats.kstest(smail.rvs(0.832555, 1, 2, size=1000),
-                        'gengamma', args=(1, 2))
+    D, p = kstest(smail(0.832555, 1, 2, size=1000), 'gengamma', args=(1, 2))
     assert p > 0.01, 'D = {}, p = {}'.format(D, p)
 
 
@@ -95,13 +54,13 @@ def test_herbel_redshift_gauss(herbel_pdf):
     cosmology = FlatLambdaCDM(H0=70, Om0=0.3, Tcmb0=2.725)
     resolution = 100
     x = np.linspace(-5, 5, resolution)
-    herbel_pdf.side_effect = [scipy.stats.norm.pdf(x)]
+    herbel_pdf.side_effect = [norm.pdf(x)]
     sample = herbel_redshift(alpha=-1.3, a_phi=-0.10268436,
                              a_m=-0.9408582, b_phi=0.00370253,
                              b_m=-20.40492365, cosmology=cosmology,
                              low=-5., high=5.0, size=1000000,
                              resolution=resolution)
-    p_value = scipy.stats.kstest(sample, 'norm')[1]
+    p_value = kstest(sample, 'norm')[1]
     assert p_value >= 0.01
 
 
@@ -126,5 +85,5 @@ def test_herbel_redshift_sampling():
         cdf = cdf / cdf[-1]
         return cdf
 
-    p_value = scipy.stats.kstest(sample, calc_cdf)[1]
+    p_value = kstest(sample, calc_cdf)[1]
     assert p_value >= 0.01