ENH: Moyal distribution

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.

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'
            ]

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)))

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:

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: