Merge pull request #353 from vhaasteren/ppf_parameters

Percentage Point Functions for parameters

enterprise/signals/parameter.py

@@ -6,6 +6,7 @@
 
 import numpy as np
 from scipy.special import erf as _erf
+import scipy.stats as sstats
 
 from enterprise.signals.selections import selection_func
 
@@ -49,9 +50,15 @@ def __init__(self, name):
             msg = "Parameter classes need to define _prior, or _logprior."
             raise AttributeError(msg)
 
+        if hasattr(self, "_ppf"):
+            self.ppf = self._ppf(name)
+        else:
+            self.ppf = None
+
         self.type = self.__class__.__name__.lower()
 
     def get_logpdf(self, value=None, **kwargs):
+        # RvH: This exception cannot be triggered
         if not isinstance(self, Parameter):
             raise TypeError("You can only call get_logpdf() on an " "instantiated (named) Parameter.")
 
@@ -66,6 +73,7 @@ def get_logpdf(self, value=None, **kwargs):
         return logpdf if self._size is None else np.sum(logpdf)
 
     def get_pdf(self, value=None, **kwargs):
+        # RvH: This exception cannot be triggered
         if not isinstance(self, Parameter):
             raise TypeError("You can only call get_pdf() on an " "instantiated (named) Parameter.")
 
@@ -87,13 +95,23 @@ def sample(self, **kwargs):
             raise AttributeError("No sampler was provided for this Parameter.")
         else:
             if self.name in kwargs:
+                # RvH: This exception cannot be triggered
                 raise ValueError("You shouldn't give me my value when you're sampling me.!")
 
             if hasattr(self, "prior"):
                 return self.prior(func=self._sampler, size=self._size, **kwargs)
             else:
                 return self.logprior(func=self._sampler, size=self._size, **kwargs)
 
+    def get_ppf(self, value=None, **kwargs):
+        if self.ppf is None:
+            raise NotImplementedError("No ppf was implemented for this Parameter.")
+
+        if value is None and "params" in kwargs:
+            value = kwargs["params"][self.name]
+
+        return self.ppf(value, **kwargs)
+
     @property
     def size(self):
         return self._size
@@ -136,7 +154,7 @@ class GPCoefficients(Parameter):
     return GPCoefficients
 
 
-def UserParameter(prior=None, logprior=None, sampler=None, size=None):
+def UserParameter(prior=None, logprior=None, sampler=None, ppf=None, size=None):
     """Class factory for UserParameter, implementing Enterprise parameters
     with arbitrary priors. The prior is specified by way of an Enterprise
     ``Function`` of the form ``prior(value, [par1, par2])``. Optionally,
@@ -147,6 +165,7 @@ def UserParameter(prior=None, logprior=None, sampler=None, size=None):
     :param prior:   parameter prior pdf, given as Enterprise ``Function``
     :param sampler: function returning a randomly sampled parameter according
                     to prior
+    :param ppf:     percentage point function (inverse cdf), for this parameter
     :param size:    length for vector parameter
     :return:        ``UserParameter`` class
     """
@@ -157,6 +176,8 @@ class UserParameter(Parameter):
             _prior = prior
         if logprior is not None:
             _logprior = logprior
+        if ppf is not None:
+            _ppf = ppf
         _sampler = None if sampler is None else staticmethod(sampler)
         _typename = "UserParameter"
 
@@ -189,6 +210,12 @@ def UniformSampler(pmin, pmax, size=None):
     return np.random.uniform(pmin, pmax, size=size)
 
 
+def UniformPPF(value, pmin, pmax):
+    """Percentage Point function for Uniform paramters."""
+
+    return sstats.uniform.ppf(value, loc=pmin, scale=pmax - pmin)
+
+
 def Uniform(pmin, pmax, size=None):
     """Class factory for Uniform parameters (with pdf(x) ~ 1/[pmax - pmin]
     inside [pmin,pmax], 0 outside. Handles vectors correctly,
@@ -204,6 +231,7 @@ def Uniform(pmin, pmax, size=None):
     class Uniform(Parameter):
         _size = size
         _prior = Function(UniformPrior, pmin=pmin, pmax=pmax)
+        _ppf = Function(UniformPPF, pmin=pmin, pmax=pmax)
         _sampler = staticmethod(UniformSampler)
         _typename = _argrepr("Uniform", pmin=pmin, pmax=pmax)
 
@@ -233,6 +261,17 @@ def NormalSampler(mu, sigma, size=None):
         return np.random.normal(mu, sigma, size=size)
 
 
+def NormalPPF(value, mu, sigma):
+    """Prior function for Normal parameters.
+    Handles scalar mu and sigma, compatible vector value/mu/sigma,
+    vector value/mu and compatible covariance matrix sigma."""
+
+    if np.ndim(sigma) >= 2:
+        raise NotImplementedError("PPF not implemented when sigma is 2D")
+
+    return sstats.norm.ppf(value, loc=mu, scale=sigma)
+
+
 def Normal(mu=0, sigma=1, size=None):
     """Class factory for Normal parameters (with pdf(x) ~ N(``mu``,``sigma``)).
     Handles vectors correctly if ``size == len(mu) == len(sigma)``,
@@ -249,6 +288,7 @@ def Normal(mu=0, sigma=1, size=None):
     class Normal(Parameter):
         _size = size
         _prior = Function(NormalPrior, mu=mu, sigma=sigma)
+        _ppf = Function(NormalPPF, mu=mu, sigma=sigma)
         _sampler = staticmethod(NormalSampler)
         _typename = _argrepr("Normal", mu=mu, sigma=sigma)
 
@@ -346,6 +386,13 @@ def LinearExpSampler(pmin, pmax, size=None):
     return np.log10(np.random.uniform(10**pmin, 10**pmax, size))
 
 
+def LinearExpPPF(value, pmin, pmax):
+    """Percentage Point function for Uniform paramters."""
+
+    ev = sstats.uniform.ppf(value, loc=10**pmin, scale=10**pmax - 10**pmin)
+    return np.log10(ev)
+
+
 def LinearExp(pmin, pmax, size=None):
     """Class factory for LinearExp parameters (with pdf(x) ~ 10^x,
     and 0 outside [``pmin``,``max``]). Handles vectors correctly
@@ -361,6 +408,7 @@ def LinearExp(pmin, pmax, size=None):
     class LinearExp(Parameter):
         _size = size
         _prior = Function(LinearExpPrior, pmin=pmin, pmax=pmax)
+        _ppf = Function(LinearExpPPF, pmin=pmin, pmax=pmax)
         _sampler = staticmethod(LinearExpSampler)
         _typename = _argrepr("LinearExp", pmin=pmin, pmax=pmax)
 
@@ -439,7 +487,6 @@ def __init__(self, name, psr=None):
 
             for kw, arg in self.func_kwargs.items():
                 if isinstance(arg, type) and issubclass(arg, (Parameter, ConstantParameter)):
-
                     # parameter name template:
                     #   pname_[signalname_][fname_]parname
                     pnames = [name, fname, kw]
@@ -582,7 +629,6 @@ def wrapper(*args, **kwargs):
                 and issubclass(arg, FunctionBase)
                 or isinstance(arg, FunctionBase)
             ):
-
                 return Function(func, **kwargs)
 
         # otherwise, we simply call the function

tests/test_parameter.py

@@ -13,10 +13,53 @@
 import numpy as np
 import scipy.stats
 
-from enterprise.signals.parameter import UniformPrior, UniformSampler, Uniform
-from enterprise.signals.parameter import NormalPrior, NormalSampler, Normal
+from enterprise.signals.parameter import Parameter, UserParameter, Function
+from enterprise.signals.parameter import UniformPrior, UniformSampler, Uniform, UniformPPF
+from enterprise.signals.parameter import NormalPrior, NormalSampler, Normal, NormalPPF
 from enterprise.signals.parameter import TruncNormalPrior, TruncNormalSampler, TruncNormal
-from enterprise.signals.parameter import LinearExpPrior, LinearExpSampler
+from enterprise.signals.parameter import LinearExpPrior, LinearExpSampler, LinearExpPPF, LinearExp
+
+
+class TestParameterExceptions(unittest.TestCase):
+    def test_missing_prior_attribute_error(self):
+        class MissingPriorParameter(Parameter):
+            pass  # Do not define _prior or _logprior
+
+        with self.assertRaises(AttributeError):
+            MissingPriorParameter("test")
+
+    def test_methods_called_on_class_type_error(self):
+        UniformClass = Uniform(pmin=0, pmax=1)
+        with self.assertRaises(TypeError):
+            UniformClass.get_logpdf()
+
+    def test_missing_sampler_attribute_error(self):
+        class MissingSamplerParameter(Parameter):
+            _prior = staticmethod(lambda x: x)
+
+            def __init__(self, name):
+                super().__init__(name)
+                self._sampler = None
+
+        missing_sampler_param = MissingSamplerParameter("test")
+        with self.assertRaises(AttributeError):
+            missing_sampler_param.sample()
+
+    def test_missing_ppf_not_implemented_error(self):
+        class MissingPPFParameter(Parameter):
+            _prior = staticmethod(lambda x: x)
+
+            def __init__(self, name):
+                super().__init__(name)
+                self.ppf = None
+
+        missing_ppf_param = MissingPPFParameter("test")
+        with self.assertRaises(NotImplementedError):
+            missing_ppf_param.get_ppf()
+
+    def test_2D_NormalPPF_error(self):
+        with self.assertRaises(NotImplementedError):
+            NormalPPF(0.0, 1.0, np.array([[1.0, 1.0], [1.0, 1.0]]))
 
 
 class TestParameter(unittest.TestCase):
@@ -35,6 +78,16 @@ def test_uniform(self):
         assert p_min < x1 < p_max, msg2
         assert type(x1) == float, msg2
 
+        msg3 = "Enterprise and scipy PPF do not match"
+        assert np.allclose(UniformPPF(x, p_min, p_max), scipy.stats.uniform.ppf(x, p_min, p_max - p_min)), msg3
+
+        # As parameter dictionary or value for Uniform instantiated object
+        unipar = Uniform(pmin=p_min, pmax=p_max)("testpar")
+        assert np.allclose(
+            unipar.get_ppf(params=dict(testpar=x)), scipy.stats.uniform.ppf(x, p_min, p_max - p_min)
+        ), msg3
+        assert np.allclose(unipar.get_ppf(x), scipy.stats.uniform.ppf(x, p_min, p_max - p_min)), msg3
+
         # vector argument
         x = np.array([0.5, 0.1])
         assert np.allclose(UniformPrior(x, p_min, p_max), scipy.stats.uniform.pdf(x, p_min, p_max - p_min)), msg1
@@ -43,9 +96,13 @@ def test_uniform(self):
         assert np.all((p_min < x1) & (x1 < p_max)), msg2
         assert x1.shape == (3,), msg2
 
+        # vector argument
+        assert np.allclose(UniformPPF(x, p_min, p_max), scipy.stats.uniform.ppf(x, p_min, p_max - p_min)), msg3
+
         # vector bounds
         p_min, p_max = np.array([0.2, 0.3]), np.array([1.1, 1.2])
         assert np.allclose(UniformPrior(x, p_min, p_max), scipy.stats.uniform.pdf(x, p_min, p_max - p_min)), msg1
+        assert np.allclose(UniformPPF(x, p_min, p_max), scipy.stats.uniform.ppf(x, p_min, p_max - p_min)), msg3
 
         x1 = UniformSampler(p_min, p_max)
         assert np.all((p_min < x1) & (x1 < p_max)), msg2
@@ -55,6 +112,33 @@ def test_uniform(self):
         assert np.all((p_min < x1) & (x1 < p_max)), msg2
         assert x1.shape == (3, 2), msg2
 
+    def test_userparameter(self):
+        """Test User-defined parameter prior, sampler, and ppf"""
+
+        # scalar
+        p_min, p_max = 0.2, 1.1
+        x = 0.5
+
+        # As parameter dictionary or value for Uniform instantiated object
+        unipar = Uniform(pmin=p_min, pmax=p_max)("testpar")
+        unipar = UserParameter(
+            prior=Function(UniformPrior, pmin=p_min, pmax=p_max),
+            sampler=staticmethod(UniformSampler),
+            ppf=Function(UniformPPF, pmin=p_min, pmax=p_max),
+        )("testpar")
+
+        msg1 = "Enterprise and scipy prior do not match"
+        assert np.allclose(
+            unipar.get_pdf(params=dict(testpar=x)), scipy.stats.uniform.pdf(x, p_min, p_max - p_min)
+        ), msg1
+        assert np.allclose(unipar.get_pdf(x), scipy.stats.uniform.pdf(x, p_min, p_max - p_min)), msg1
+
+        msg2 = "Enterprise and scipy PPF do not match"
+        assert np.allclose(
+            unipar.get_ppf(params=dict(testpar=x)), scipy.stats.uniform.ppf(x, p_min, p_max - p_min)
+        ), msg2
+        assert np.allclose(unipar.get_ppf(x), scipy.stats.uniform.ppf(x, p_min, p_max - p_min)), msg2
+
     def test_linearexp(self):
         """Test LinearExp parameter prior and sampler."""
 
@@ -68,6 +152,19 @@ def test_linearexp(self):
         msg1b = "Scalar sampler out of range"
         assert p_min <= x <= p_max, msg1b
 
+        msg1c = "Scalar PPF does not match"
+        x = 0.5
+        assert np.allclose(
+            LinearExpPPF(x, p_min, p_max), np.log10(10**p_min + x * (10**p_max - 10**p_min))
+        ), msg1c
+
+        # As parameter dictionary or value for Uniform instantiated object
+        lepar = LinearExp(pmin=p_min, pmax=p_max)("testpar")
+        assert np.allclose(
+            lepar.get_ppf(params=dict(testpar=x)), np.log10(10**p_min + x * (10**p_max - 10**p_min))
+        ), msg1
+        assert np.allclose(lepar.get_ppf(x), np.log10(10**p_min + x * (10**p_max - 10**p_min))), msg1
+
         # vector argument
         x = np.array([0, 1.5, 2.5])
         msg2 = "Vector-argument prior does not match"
@@ -79,6 +176,12 @@ def test_linearexp(self):
         msg2b = "Vector-argument sampler out of range"
         assert np.all((p_min < x) & (x < p_max)), msg2b
 
+        x = np.array([0.5, 0.75])
+        msg2c = "Vector-argument PPF does not match"
+        assert np.allclose(
+            LinearExpPPF(x, p_min, p_max), np.log10(10**p_min + x * (10**p_max - 10**p_min))
+        ), msg2c
+
         # vector bounds
         p_min, p_max = np.array([0, 1]), np.array([2, 3])
         x = np.array([1, 2])
@@ -88,6 +191,14 @@ def test_linearexp(self):
             np.array([10**1 / (10**2 - 10**0), 10**2 / (10**3 - 10**1)]) * np.log(10),
         ), msg3
 
+        # Vector PPF
+        x = np.array([0.5, 0.75])
+        p_min, p_max = np.array([0, 1]), np.array([2, 3])
+        msg3c = "Vector-argument PPF+bounds does not match"
+        assert np.allclose(
+            LinearExpPPF(x, p_min, p_max), np.log10(10**p_min + x * (10**p_max - 10**p_min))
+        ), msg3c
+
     def test_normal(self):
         """Test Normal parameter prior and sampler for various combinations of scalar and vector arguments."""
 
@@ -105,6 +216,9 @@ def test_normal(self):
         # this should almost never fail
         assert -5 < (x1 - mu) / sigma < 5, msg2
 
+        msg3 = "Enterprise and scipy PPF do not match"
+        assert np.allclose(NormalPPF(x, mu, sigma), scipy.stats.norm.ppf(x, loc=mu, scale=sigma)), msg3
+
         # vector argument
         x = np.array([-0.2, 0.1, 0.5])
 
@@ -118,6 +232,9 @@ def test_normal(self):
         )
         assert x1.shape == x2.shape, msg2
 
+        x = np.array([0.1, 0.25, 0.65])
+        assert np.allclose(NormalPPF(x, mu, sigma), scipy.stats.norm.ppf(x, loc=mu, scale=sigma)), msg3
+
         # vector bounds; note the different semantics from `NormalPrior`,
         # which returns a vector consistently with `UniformPrior`
         mu, sigma = np.array([0.1, 0.15, 0.2]), np.array([2, 1, 2])
@@ -199,4 +316,4 @@ def test_metaparam(self):
 
         paramA = TruncNormal(mu, sigma, pmin, pmax)("A")
         xs = np.array([-3.5, 3.5])
-        assert np.alltrue(paramA.get_pdf(xs, mu=mu.sample()) == zeros), msg4
+        assert np.all(paramA.get_pdf(xs, mu=mu.sample()) == zeros), msg4