diff --git a/enterprise/signals/parameter.py b/enterprise/signals/parameter.py index 38e4a4e2..34c8f3be 100644 --- a/enterprise/signals/parameter.py +++ b/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,6 +95,7 @@ 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"): @@ -94,6 +103,15 @@ def sample(self, **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 diff --git a/tests/test_parameter.py b/tests/test_parameter.py index b706b834..c549ea27 100644 --- a/tests/test_parameter.py +++ b/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