From 8abcee8da7fe0c74812b3115181904d1ca6f18b9 Mon Sep 17 00:00:00 2001 From: Fearghus Robert Keeble Date: Mon, 6 Apr 2020 16:29:05 +0100 Subject: [PATCH 1/3] Added Moyal distribution to continuous class. Added tests for Moyal distribution. --- pymc3/distributions/__init__.py | 2 + pymc3/distributions/continuous.py | 134 ++++++++++++++++++++++- pymc3/tests/test_distributions.py | 8 +- pymc3/tests/test_distributions_random.py | 11 ++ 4 files changed, 152 insertions(+), 3 deletions(-) diff --git a/pymc3/distributions/__init__.py b/pymc3/distributions/__init__.py index a9e69fbacdc..bfecad95ef9 100644 --- a/pymc3/distributions/__init__.py +++ b/pymc3/distributions/__init__.py @@ -48,6 +48,7 @@ from .continuous import LogitNormal from .continuous import Interpolated from .continuous import Rice +from .continuous import Moyal from .discrete import Binomial from .discrete import BetaBinomial @@ -170,6 +171,7 @@ 'Interpolated', 'Bound', 'Rice', + 'Moyal', 'Simulator', 'fast_sample_posterior_predictive' ] diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py index 55befc0201e..51d3bc68bfd 100644 --- a/pymc3/distributions/continuous.py +++ b/pymc3/distributions/continuous.py @@ -40,7 +40,7 @@ 'HalfCauchy', 'Gamma', 'Weibull', 'HalfStudentT', 'Lognormal', 'ChiSquared', 'HalfNormal', 'Wald', 'Pareto', 'InverseGamma', 'ExGaussian', 'VonMises', 'SkewNormal', 'Triangular', 'Gumbel', - 'Logistic', 'LogitNormal', 'Interpolated', 'Rice'] + 'Logistic', 'LogitNormal', 'Interpolated', 'Rice', 'Moyal'] class PositiveContinuous(Continuous): @@ -1946,7 +1946,7 @@ class StudentT(Continuous): plt.show() ======== ======================== - Support :math:``x \in \mathbb{R}`` + Support :math:`x \in \mathbb{R}` ======== ======================== Parameters @@ -4381,3 +4381,133 @@ def logp(self, value): TensorVariable """ return tt.log(self.interp_op(value) / self.Z) + + +class Moyal(Continuous): + R""" + Moyal log-likelihood. + + The pdf of this distribution is + + .. math:: + + f(x \mid \mu,\sigma) = \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{1}{2}\left(z + e^{-z}\right)}, + + where + + .. math:: + + z = \frac{x-\mu}{\sigma}. + + .. plot:: + + import matplotlib.pyplot as plt + import numpy as np + import scipy.stats as st + plt.style.use('seaborn-darkgrid') + x = np.linspace(-10, 20, 200) + mus = [-1., 0., 4.] + sigmas = [2., 2., 4.] + for mu, sigma in zip(mus, sigmas): + pdf = st.moyal.pdf(x, loc=mu, scale=sigma) + plt.plot(x, pdf, label=r'$\mu$ = {}, $\sigma$ = {}'.format(mu, sigma)) + plt.xlabel('x', fontsize=12) + plt.ylabel('f(x)', fontsize=12) + plt.legend(loc=1) + plt.show() + + ======== ============================================================== + Support :math:`x \in (-\infty, \infty)` + Mean :math:`\mu + \sigma\left(\gamma + \log 2\right)`, where :math:`\gamma` is the Euler-Mascheroni constant + Variance :math:`\frac{\pi^{2}}{2}\sigma^{2}` + ======== ============================================================== + + Parameters + ---------- + mu: float + Location parameter. + sigma: float + Scale parameter (sigma > 0). + """ + + def __init__(self, mu=0, sigma=1., *args, **kwargs): + self.mu = tt.as_tensor_variable(floatX(mu)) + self.sigma = tt.as_tensor_variable(floatX(sigma)) + + assert_negative_support(sigma, 'sigma', 'Moyal') + + self.mean = self.mu + self.sigma * (np.euler_gamma + tt.log(2)) + self.median = self.mu - self.sigma * tt.log(2 * tt.erfcinv(1 / 2)**2) + self.mode = self.mu + self.variance = (np.pi**2 / 2.0) * self.sigma**2 + + super().__init__(*args, **kwargs) + + def random(self, point=None, size=None): + """ + Draw random values from Moyal distribution. + + Parameters + ---------- + point: dict, optional + Dict of variable values on which random values are to be + conditioned (uses default point if not specified). + size: int, optional + Desired size of random sample (returns one sample if not + specified). + + Returns + ------- + array + """ + mu, sigma = draw_values([self.mu, self.sigma], point=point, size=size) + return generate_samples(stats.moyal.rvs, loc=mu, scale=sigma, + dist_shape=self.shape, + size=size) + + def logp(self, value): + """ + Calculate log-probability of Moyal distribution at specified value. + + Parameters + ---------- + value: numeric + Value(s) for which log-probability is calculated. If the log probabilities for multiple + values are desired the values must be provided in a numpy array or theano tensor + + Returns + ------- + TensorVariable + """ + scaled = (value - self.mu) / self.sigma + return bound((-(1 / 2) * (scaled + tt.exp(-scaled)) + - tt.log(self.sigma) + - (1 / 2) * tt.log(2 * np.pi)), self.sigma > 0) + + def _repr_latex_(self, name=None, dist=None): + if dist is None: + dist = self + sigma = dist.sigma + mu = dist.mu + name = r'\text{%s}' % name + return r'${} \sim \text{{Moyal}}(\mathit{{mu}}={},~\mathit{{sigma}}={})$'.format(name, + get_variable_name(mu), + get_variable_name(sigma)) + + def logcdf(self, value): + """ + Compute the log of the cumulative distribution function for Moyal distribution + at the specified value. + + Parameters + ---------- + value: numeric + Value(s) for which log CDF is calculated. If the log CDF for multiple + values are desired the values must be provided in a numpy array or theano tensor. + + Returns + ------- + TensorVariable + """ + scaled = (value - self.mu) / self.sigma + return tt.log(tt.erfc(tt.exp(-scaled / 2) * (2**-0.5))) diff --git a/pymc3/tests/test_distributions.py b/pymc3/tests/test_distributions.py index 93534054515..65138ce62ba 100644 --- a/pymc3/tests/test_distributions.py +++ b/pymc3/tests/test_distributions.py @@ -30,7 +30,7 @@ Bound, Uniform, Triangular, Binomial, SkewNormal, DiscreteWeibull, Gumbel, Logistic, OrderedLogistic, LogitNormal, Interpolated, ZeroInflatedBinomial, HalfFlat, AR1, KroneckerNormal, Rice, - Kumaraswamy + Kumaraswamy, Moyal ) from ..distributions import continuous @@ -1213,6 +1213,12 @@ def test_rice(self): self.pymc3_matches_scipy(Rice, Rplus, {'b': Rplus, 'sigma': Rplusbig}, lambda value, b, sigma: sp.rice.logpdf(value, b=b, loc=0, scale=sigma)) + def test_moyal(self): + self.pymc3_matches_scipy(Moyal, R, {'mu': R, 'sigma': Rplusbig}, + lambda value, mu, sigma: sp.moyal.logpdf(value, mu, sigma)) + self.check_logcdf(Moyal, R, {'mu': R, 'sigma': Rplusbig}, + lambda value, mu, sigma: sp.moyal.logcdf(value, mu, sigma)) + @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32") def test_interpolated(self): for mu in R.vals: diff --git a/pymc3/tests/test_distributions_random.py b/pymc3/tests/test_distributions_random.py index 393af88ad83..5da7c4ce5ee 100644 --- a/pymc3/tests/test_distributions_random.py +++ b/pymc3/tests/test_distributions_random.py @@ -466,7 +466,12 @@ class TestGeometric(BaseTestCases.BaseTestCase): distribution = pm.Geometric params = {'p': 0.5} + +class TestMoyal(BaseTestCases.BaseTestCase): + distribution = pm.Moyal + params = {'mu': 0., 'sigma': 1.} + class TestCategorical(BaseTestCases.BaseTestCase): distribution = pm.Categorical params = {'p': np.ones(BaseTestCases.BaseTestCase.shape)} @@ -821,6 +826,12 @@ def ref_rand(size, mu, sigma): return expit(st.norm.rvs(loc=mu, scale=sigma, size=size)) pymc3_random(pm.LogitNormal, {'mu': R, 'sigma': Rplus}, ref_rand=ref_rand) + def test_moyal(self): + def ref_rand(size, mu, sigma): + return st.moyal.rvs(loc=mu, scale=sigma, size=size) + pymc3_random(pm.Moyal, {'mu': R, 'sigma': Rplus}, ref_rand=ref_rand) + + @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32") def test_interpolated(self): for mu in R.vals: From 50a85f8925ac820828e434485cd6310d971fb1b0 Mon Sep 17 00:00:00 2001 From: Fearghus Robert Keeble Date: Mon, 6 Apr 2020 16:55:13 +0100 Subject: [PATCH 2/3] Release notes updated. --- RELEASE-NOTES.md | 1 + 1 file changed, 1 insertion(+) diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md index a4f4a950999..c7c80375474 100644 --- a/RELEASE-NOTES.md +++ b/RELEASE-NOTES.md @@ -11,6 +11,7 @@ - `sample_posterior_predictive` can now feed on `xarray.Dataset` - e.g. from `InferenceData.posterior`. (see [#3846](https://github.com/pymc-devs/pymc3/pull/3846)) - `SamplerReport` (`MultiTrace.report`) now has properties `n_tune`, `n_draws`, `t_sampling` for increased convenience (see [#3827](https://github.com/pymc-devs/pymc3/pull/3827)) - `pm.sample` now has support for adapting dense mass matrix using `QuadPotentialFullAdapt` (see [#3596](https://github.com/pymc-devs/pymc3/pull/3596), [#3705](https://github.com/pymc-devs/pymc3/pull/3705) and [#3858](https://github.com/pymc-devs/pymc3/pull/3858)) +- `Moyal` distribution added (see [#3870](https://github.com/pymc-devs/pymc3/pull/3870)). ### Maintenance - Remove `sample_ppc` and `sample_ppc_w` that were deprecated in 3.6. From d3daef9bcb0f81900db164b0b0c75e16801eddf9 Mon Sep 17 00:00:00 2001 From: Fearghus Robert Keeble Date: Tue, 7 Apr 2020 10:14:42 +0100 Subject: [PATCH 3/3] Added float32 exception for Moyal dist in tests. --- pymc3/tests/test_distributions.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pymc3/tests/test_distributions.py b/pymc3/tests/test_distributions.py index 65138ce62ba..d8586d421c5 100644 --- a/pymc3/tests/test_distributions.py +++ b/pymc3/tests/test_distributions.py @@ -1213,11 +1213,12 @@ def test_rice(self): self.pymc3_matches_scipy(Rice, Rplus, {'b': Rplus, 'sigma': Rplusbig}, lambda value, b, sigma: sp.rice.logpdf(value, b=b, loc=0, scale=sigma)) + @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32") def test_moyal(self): self.pymc3_matches_scipy(Moyal, R, {'mu': R, 'sigma': Rplusbig}, - lambda value, mu, sigma: sp.moyal.logpdf(value, mu, sigma)) + lambda value, mu, sigma: floatX(sp.moyal.logpdf(value, mu, sigma))) self.check_logcdf(Moyal, R, {'mu': R, 'sigma': Rplusbig}, - lambda value, mu, sigma: sp.moyal.logcdf(value, mu, sigma)) + lambda value, mu, sigma: floatX(sp.moyal.logcdf(value, mu, sigma))) @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32") def test_interpolated(self):