Repository: skypyproject/skypy
Pull request #211 (2020-07-09 13:11:10)

Sampling of the galaxy stellar masses
## Description
This module samples the galaxy stellar masses from a Schechter function according to the model of [1].
## Checklist
- [ ] Follow the [Contributor Guidelines](https://github.com/skypyproject/skypy/blob/master/CONTRIBUTING.md)
- [ ] Write unit tests
- [ ] Write documentation strings
- [ ] Assign someone from your working team to review this pull request
- [ ] Assign someone from the infrastructure team to review this pull request

## References
[1] Peng et al. (2010)

philipp128: > There is no redshift dependence in this model so it doesn't make sense for this to be a function of redshift and the failing test is meaningless. However, we should generalise the function to accept m_star as an array. Then it would make sense for the function to raise an exception if size != m_star.size. I've given suggestions for how to make this change but you will need to update the rest of your examples and tests as well to make them consistent.

True. I am not sure whether it makes sense to make `m_star` array-like but not `alpha`?
rrjbca: This suggestion is equivalent to what we did for the luminosity function where L* was an array and alpha was constant. The only difference in that case was that instead of passing L* to the function directly we passed an array of redshifts and calculated L* at each redshift inside the function.
skypybot: Hello @philipp128! Thanks for updating 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_stellar_mass.py`](https://github.com/skypyproject/skypy/blob/3933708666dc7d8ef146b21b873d7a9c06b61df8/skypy/galaxy/tests/test_stellar_mass.py):

> [Line 53:80](https://github.com/skypyproject/skypy/blob/3933708666dc7d8ef146b21b873d7a9c06b61df8/skypy/galaxy/tests/test_stellar_mass.py#L53): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (83 > 79 characters)


rrjbca: We also need a test where `size = None` and m_star is an array to improve the test coverage. E.g. something that checks `m_star.shape == schechter_mass(..., m_star, ...).shape`
ntessore: Looks okay now, but the new module is not included in the documentation yet.
philipp128: > Looks okay now, but the new module is not included in the documentation yet.

How do you do this again? Which file is it?

Patch:
diff --git a/skypy/galaxy/stellar_mass.py b/skypy/galaxy/stellar_mass.py
new file mode 100644
index 0000000..1efc83a
--- /dev/null
+++ b/skypy/galaxy/stellar_mass.py
@@ -0,0 +1,75 @@
+r"""Models of galaxy stellar mass.
+
+"""
+
+import numpy as np
+
+from skypy.utils.random import schechter
+
+
+__all__ = [
+    'schechter_smf',
+]
+
+
+def schechter_smf(alpha, m_star, x_min, x_max, resolution=100, size=None):
+    r""" Stellar masses following the Schechter mass function [1]_.
+
+    Parameters
+    ----------
+    alpha : float
+        The alpha parameter in the Schechter stellar mass function.
+    m_star : array-like
+        Characteristic stellar mass M_*.
+    size: int, optional
+         Output shape of stellar mass samples. If size is None and m_star
+         is a scalar, a single sample is returned. If size is None and
+         m_star is an array, an array of samples is returned with the same
+         shape as m_star.
+    x_min, x_max : float
+        Lower and upper bounds for the random variable x in units of M_*.
+    resolution : int, optional
+        Resolution of the inverse transform sampling spline. Default is 100.
+
+    Returns
+    -------
+    stellar mass : array_like
+        Drawn stellar masses from the Schechter stellar mass function in units
+        of the solar mass
+
+    Notes
+    -----
+     The stellar mass probability distribution (pdf) follows a Schechter
+    profile of the form
+
+    .. math::
+        \Phi(M) = \frac{1}{M_*} \left(\frac{M}{M_*}\right)^\alpha
+            \exp\left(-\frac{M}{M_*}\right) \;.
+
+    From this pdf one can sample the stellar masses.
+
+
+    References
+    ----------
+    .. [1] Mo, H., Van den Bosch, F., & White, S. (2010). Galaxy Formation and
+    Evolution. Cambridge: Cambridge University Press.
+    doi:10.1017/CBO9780511807244
+
+    Examples
+    --------
+    >>> from skypy.galaxy import stellar_mass
+
+    Sample 100 stellar masses values at redshift z = 1.0 with alpha = -1.4,
+    m_star = 10**10.67, x_min = 0.0002138 and x_max = 213.8
+
+    >>> masses = stellar_mass.schechter_smf(-1.4, 10**10.67, 0.0002138,
+    ...                               213.8, size=100)
+
+    """
+
+    if size is None and np.shape(m_star):
+        size = np.shape(m_star)
+
+    x_sample = schechter(alpha, x_min, x_max, resolution=resolution, size=size)
+
+    return m_star * x_sample


Test patch:
diff --git a/skypy/galaxy/tests/test_stellar_mass.py b/skypy/galaxy/tests/test_stellar_mass.py
new file mode 100644
index 0000000..1730872
--- /dev/null
+++ b/skypy/galaxy/tests/test_stellar_mass.py
@@ -0,0 +1,62 @@
+import numpy as np
+import scipy.stats
+import scipy.integrate
+from scipy.special import gammaln
+import pytest
+
+from skypy.galaxy import stellar_mass
+from skypy.utils import special
+
+
+def test_exponential_distribution():
+    # When alpha=0, M*=1 and x_min~0 we get a truncated exponential
+    q_max = 1e2
+    sample = stellar_mass.schechter_smf(0, 1, size=1000,
+                                        x_min=1e-10, x_max=q_max,
+                                        resolution=1000)
+    d, p_value = scipy.stats.kstest(sample, 'truncexpon', args=(q_max,))
+    assert p_value >= 0.01
+
+
+def test_stellar_masses():
+    # Test that error is returned if m_star input is an array but size !=
+    # None and size != m_star,size
+    with pytest.raises(ValueError):
+        stellar_mass.schechter_smf(-1.4, np.array([1e10, 2e10]), 2, 1, size=3)
+
+    # Test that an array with the sme shape as m_star is returned if m_star is
+    # an array and size = None
+    m_star = np.array([1e10, 2e10])
+    sample = stellar_mass.schechter_smf(-1.4, m_star, 2, 1, size=None)
+    assert m_star.shape == sample.shape
+
+    # Test that sampling corresponds to sampling from the right pdf.
+    # For this, we sample an array of luminosities for redshift z = 1.0 and we
+    # compare it to the corresponding cdf.
+
+    def calc_pdf(m, alpha, mass_star, mass_min, mass_max):
+        lg = gammaln(alpha + 1)
+        c = np.fabs(special.gammaincc(alpha + 1, mass_min / mass_star))
+        d = np.fabs(special.gammaincc(alpha + 1, mass_max / mass_star))
+        norm = np.exp(lg) * (c - d)
+        return 1. / mass_star * np.power(m / mass_star, alpha) * \
+            np.exp(-m / mass_star) / norm
+
+    def calc_cdf(m):
+        alpha = -1.4
+        mass_star = 10 ** 10.67
+        mass_min = 10 ** 7
+        mass_max = 10 ** 13
+        pdf = calc_pdf(m, alpha, mass_star, mass_min, mass_max)
+        cdf = scipy.integrate.cumtrapz(pdf, m, initial=0)
+        cdf = cdf / cdf[-1]
+        return cdf
+
+    m_star = 10 ** 10.67
+    m_min = 10 ** 7
+    m_max = 10 ** 13
+    sample = stellar_mass.schechter_smf(-1.4, m_star, m_min / m_star,
+                                        m_max / m_star, size=1000,
+                                        resolution=100)
+    p_value = scipy.stats.kstest(sample, calc_cdf)[1]
+    assert p_value >= 0.01