diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md index a4f4a95099..c7c8037547 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. diff --git a/pymc3/distributions/__init__.py b/pymc3/distributions/__init__.py index a9e69fbacd..bfecad95ef 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 55befc0201..51d3bc68bf 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 9353405451..d8586d421c 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,13 @@ 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: floatX(sp.moyal.logpdf(value, mu, sigma))) + self.check_logcdf(Moyal, R, {'mu': R, 'sigma': Rplusbig}, + 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): for mu in R.vals: diff --git a/pymc3/tests/test_distributions_random.py b/pymc3/tests/test_distributions_random.py index 393af88ad8..5da7c4ce5e 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: