Repository: skypyproject/skypy
Pull request #199 (2020-06-23 13:08:09)

simplify beta ellipticity generation
## Description

Simplifies the beta ellipticity model to only generate random variates using `numpy.random.beta`. If that's undesirable please close. Part of the #198 fix.

## 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

skypybot: Hello @ntessore! Thanks for opening this PR. We checked the lines you've touched for [PEP 8](https://www.python.org/dev/peps/pep-0008) issues, and found:



* In the file [`skypy/galaxy/tests/test_ellipticity.py`](https://github.com/skypyproject/skypy/blob/f80e2bf07cd5017ca8b5cbdf1c384d8f6cb1d968/skypy/galaxy/tests/test_ellipticity.py):

> [Line 10:80](https://github.com/skypyproject/skypy/blob/f80e2bf07cd5017ca8b5cbdf1c384d8f6cb1d968/skypy/galaxy/tests/test_ellipticity.py#L10): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (80 > 79 characters)


codecov[bot]: # [Codecov](https://codecov.io/gh/skypyproject/skypy/pull/199?src=pr&el=h1) Report
> Merging [#199](https://codecov.io/gh/skypyproject/skypy/pull/199?src=pr&el=desc) into [master](https://codecov.io/gh/skypyproject/skypy/commit/573a0a0239538a64efd8194d5e3944cca4b84253&el=desc) will **increase** coverage by `0.18%`.
> The diff coverage is `100.00%`.

[![Impacted file tree graph](https://codecov.io/gh/skypyproject/skypy/pull/199/graphs/tree.svg?width=650&height=150&src=pr&token=9jNO45edcU)](https://codecov.io/gh/skypyproject/skypy/pull/199?src=pr&el=tree)

```diff
@@            Coverage Diff             @@
##           master     #199      +/-   ##
==========================================
+ Coverage   99.61%   99.80%   +0.18%     
==========================================
  Files          19       19              
  Lines         517      500      -17     
==========================================
- Hits          515      499      -16     
+ Misses          2        1       -1     
```


| [Impacted Files](https://codecov.io/gh/skypyproject/skypy/pull/199?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [skypy/galaxy/ellipticity.py](https://codecov.io/gh/skypyproject/skypy/pull/199/diff?src=pr&el=tree#diff-c2t5cHkvZ2FsYXh5L2VsbGlwdGljaXR5LnB5) | `100.00% <100.00%> (+2.43%)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/skypyproject/skypy/pull/199?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/skypyproject/skypy/pull/199?src=pr&el=footer). Last update [573a0a0...f80e2bf](https://codecov.io/gh/skypyproject/skypy/pull/199?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

ntessore: @rrjbca It would probably be nicer to have realistic parameters in the example to get a realistic [plot in the docs](https://996-221432358-gh.circle-artifacts.com/0/html/api/skypy.galaxy.ellipticity.beta_ellipticity.html#skypy.galaxy.ellipticity.beta_ellipticity). I cannot seem to find any, though?

Patch:
diff --git a/skypy/galaxy/ellipticity.py b/skypy/galaxy/ellipticity.py
index 7a7a732..6f0eb92 100644
--- a/skypy/galaxy/ellipticity.py
+++ b/skypy/galaxy/ellipticity.py
@@ -13,16 +13,7 @@ __all__ = [
 ]
 
 
-def _ellipticity_parameterization(beta_method):
-    def reparameterize(self, *args):
-        *beta_args, e_ratio, e_sum = args
-        a = e_sum * e_ratio
-        b = e_sum * (1.0 - e_ratio)
-        return beta_method(self, *beta_args, a, b)
-    return reparameterize
-
-
-class beta_ellipticity_gen(stats._continuous_distns.beta_gen):
+def beta_ellipticity(e_ratio, e_sum, size=None):
     r'''Galaxy ellipticities sampled from a reparameterized beta distribution.
 
     The ellipticities follow a beta distribution parameterized by
@@ -57,46 +48,19 @@ class beta_ellipticity_gen(stats._continuous_distns.beta_gen):
     --------
     >>> from skypy.galaxy.ellipticity import beta_ellipticity
 
-    Sample 10 random variates from the Kacprzak model with
+    Sample 10000 random variates from the Kacprzak model with
     :math:`e_{\rm ratio} = 0.5` and :math:`e_{\rm sum} = 1.0`:
 
-    >>> ellipticity = beta_ellipticity.rvs(0.5, 1.0, size=10)
-
-    Fix distribution parameters for repeated use:
+    >>> ellipticity = beta_ellipticity(0.5, 1.0, size=10000)
 
-    >>> ellipticity_distribution = beta_ellipticity(0.5, 1.0)
-    >>> ellipticity_distribution.mean()
-    0.5
-    >>> ellipticity = ellipticity_distribution.rvs(size=10)
     '''
 
-    @_ellipticity_parameterization
-    def _rvs(self, *args):
-        return super()._rvs(*args)
-
-    @_ellipticity_parameterization
-    def _logpdf(self, *args):
-        return super()._logpdf(*args)
-
-    @_ellipticity_parameterization
-    def _cdf(self, *args):
-        return super()._cdf(*args)
-
-    @_ellipticity_parameterization
-    def _ppf(self, *args):
-        return super()._ppf(*args)
-
-    @_ellipticity_parameterization
-    def _stats(self, *args):
-        return super()._stats(*args)
-
-    def fit(self, data, *args, **kwargs):
-        return super(stats._continuous_distns.beta_gen, self).fit(
-            data, *args, **kwargs)
-
+    # convert to beta distribution parameters
+    a = e_sum * e_ratio
+    b = e_sum * (1.0 - e_ratio)
 
-beta_ellipticity = beta_ellipticity_gen(a=0.0, b=1.0, name='beta_ellipticity',
-                                        shapes="e_ratio e_sum")
+    # sample from the beta distribution
+    return np.random.beta(a, b, size)
 
 
 def ryden04(mu_gamma, sigma_gamma, mu, sigma, size=None):


Test patch:
diff --git a/skypy/galaxy/tests/test_ellipticity.py b/skypy/galaxy/tests/test_ellipticity.py
index 605a420..d077da5 100644
--- a/skypy/galaxy/tests/test_ellipticity.py
+++ b/skypy/galaxy/tests/test_ellipticity.py
@@ -1,80 +1,39 @@
-from astropy.units import allclose, isclose
 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)
 
 
 def test_beta_ellipticity():
 
     from skypy.galaxy.ellipticity import beta_ellipticity
 
-    # Initialise a randomised ellipticity distribution, an equivalent beta
-    # distribution and special cases where the ellipticity distribution is
-    # equivalent to a uniform distribution and an arcsine distribution.
+    # randomised ellipticity distribution with beta distribution parameters a,b
+    # and the equivalent reparametrisation
     a, b = np.random.lognormal(size=2)
-    args = (a / (a + b), a + b)
-    beta_dist = stats.beta(a, b)
-    ellipticity_dist = beta_ellipticity(*args)
-    ellipticity_uniform = beta_ellipticity(0.5, 2.0)
-    ellipticity_arcsine = beta_ellipticity(0.5, 1.0)
-
-    # Range of input values spanning the support of the distributions
-    x = np.linspace(0, 1, 100)
-
-    # Check basic properties of distribution implementation
-    check_normalization(beta_ellipticity, args, 'beta_ellipticity')
-    check_edge_support(beta_ellipticity, args)
-    check_random_state_property(beta_ellipticity, args)
-    check_pickling(beta_ellipticity, args)
-
-    # Check distribution moments
-    m, v, s, k = ellipticity_dist.stats(moments='mvsk')
-    check_mean_expect(beta_ellipticity, args, m, 'beta_ellipticity')
-    check_var_expect(beta_ellipticity, args, m, v, 'beta_ellipticity')
-    check_skew_expect(beta_ellipticity, args, m, v, s, 'beta_ellipticity')
-    check_kurt_expect(beta_ellipticity, args, m, v, k, 'beta_ellipticity')
-    check_moment(beta_ellipticity, args, m, v, 'beta_ellipticity')
-
-    # Compare ellipticity distribution functions (e.g. pdf, cdf...) against
-    # functions for an equivalent beta distribution
-    assert allclose(ellipticity_dist.pdf(x), beta_dist.pdf(x))
-    assert allclose(ellipticity_dist.logpdf(x), beta_dist.logpdf(x))
-    assert allclose(ellipticity_dist.cdf(x), beta_dist.cdf(x))
-    assert allclose(ellipticity_dist.logcdf(x), beta_dist.logcdf(x))
-    assert allclose(ellipticity_dist.ppf(x), beta_dist.ppf(x))
-    assert allclose(ellipticity_dist.sf(x), beta_dist.sf(x))
-    assert allclose(ellipticity_dist.logsf(x), beta_dist.logsf(x))
-    assert allclose(ellipticity_dist.isf(x), beta_dist.isf(x))
-    assert isclose(ellipticity_dist.entropy(), beta_dist.entropy())
-    assert isclose(ellipticity_dist.median(), beta_dist.median())
-    assert isclose(ellipticity_dist.std(), beta_dist.std())
-    assert allclose(ellipticity_dist.interval(x), beta_dist.interval(x))
+    e_ratio, e_sum = (a / (a + b), a + b)
 
     # Test scalar output
-    assert np.isscalar(ellipticity_dist.rvs())
+    assert np.isscalar(beta_ellipticity(e_ratio, e_sum))
 
     # Test array output
-    assert ellipticity_dist.rvs(size=10).shape == (10,)
+    assert beta_ellipticity(e_ratio, e_sum, size=10).shape == (10,)
 
     # Test broadcast output
-    e_ratio = 0.5 * np.ones((13, 1, 5))
-    e_sum = 0.5 * np.ones((7, 5))
-    rvs = beta_ellipticity.rvs(e_ratio, e_sum)
-    assert rvs.shape == np.broadcast(e_ratio, e_sum).shape
+    e_ratio2 = 0.5 * np.ones((13, 1, 5))
+    e_sum2 = 0.5 * np.ones((7, 5))
+    rvs = beta_ellipticity(e_ratio2, e_sum2)
+    assert rvs.shape == np.broadcast(e_ratio2, e_sum2).shape
 
     # Kolmogorov-Smirnov test comparing ellipticity and beta distributions
-    D, p = stats.kstest(ellipticity_dist.rvs, beta_dist.cdf, N=1000)
+    D, p = stats.kstest(beta_ellipticity(e_ratio, e_sum, size=1000), 'beta',
+                        args=(a, b))
     assert p > 0.01, 'D = {}, p = {}'.format(D, p)
 
     # Kolmogorov-Smirnov test comparing ellipticity and uniform distributions
-    D, p = stats.kstest(ellipticity_uniform.rvs, 'uniform', N=1000)
+    D, p = stats.kstest(beta_ellipticity(0.5, 2.0, size=1000), 'uniform')
     assert p > 0.01, 'D = {}, p = {}'.format(D, p)
 
     # Kolmogorov-Smirnov test comparing ellipticity and arcsine distributions
-    D, p = stats.kstest(ellipticity_arcsine.rvs, 'arcsine', N=1000)
+    D, p = stats.kstest(beta_ellipticity(0.5, 1.0, size=1000), 'arcsine')
     assert p > 0.01, 'D = {}, p = {}'.format(D, p)