Add logistic distribution (#2588)

Author: AustinRochford

pymc3/distributions/__init__.py

@@ -26,6 +26,7 @@
 from .continuous import SkewNormal
 from .continuous import Triangular
 from .continuous import Gumbel
+from .continuous import Logistic
 from .continuous import Interpolated
 
 from .discrete import Binomial
@@ -141,6 +142,7 @@
            'Triangular',
            'DiscreteWeibull',
            'Gumbel',
+           'Logistic',
            'Interpolated',
            'Bound',
            ]

pymc3/distributions/continuous.py

@@ -28,7 +28,7 @@
            'Laplace', 'StudentT', 'Cauchy', 'HalfCauchy', 'Gamma', 'Weibull',
            'HalfStudentT', 'StudentTpos', 'Lognormal', 'ChiSquared',
            'HalfNormal', 'Wald', 'Pareto', 'InverseGamma', 'ExGaussian',
-           'VonMises', 'SkewNormal', 'Interpolated']
+           'VonMises', 'SkewNormal', 'Logistic', 'Interpolated']
 
 
 class PositiveContinuous(Continuous):
@@ -1931,6 +1931,83 @@ def _repr_latex_(self, name=None, dist=None):
                                                                 get_variable_name(beta))
 
 
+class Logistic(Continuous):
+    R"""
+    Logistic log-likelihood.
+
+    .. math::
+
+       f(x \mid \mu, s) =
+           \frac{\exp\left(-\frac{x - \mu}{s}\right)}{s \left(1 + \exp\left(-\frac{x - \mu}{s}\right)\right)^2}
+
+    ========  ==========================================
+    Support   :math:`x \in \mathbb{R}`
+    Mean      :math:`\mu`
+    Variance  :math:`\frac{s^2 \pi^2}{3}`
+    ========  ==========================================
+
+    .. plot::
+
+        import matplotlib.pyplot as plt
+        import numpy as np
+        import scipy.stats as st
+        x = np.linspace(-5.0, 5.0, 1000)
+        fig, ax = plt.subplots()
+        f = lambda mu, s : st.logistic.pdf(x, loc=mu, scale=s)
+        plot_pdf = lambda a, b : ax.plot(x, f(a,b), label=r'$\mu$={0}, $s$={1}'.format(a,b))
+        plot_pdf(0.0, 0.4)
+        plot_pdf(0.0, 1.0)
+        plot_pdf(0.0, 2.0)
+        plot_pdf(-2.0, 0.4)
+        plt.legend(loc='upper right', frameon=False)
+        ax.set(xlim=[-5,5], ylim=[0,1.2], xlabel='x', ylabel='f(x)')
+        plt.show()
+
+    Parameters
+    ----------
+    mu : float
+        Mean.
+    s : float
+        Scale (s > 0).
+    """
+    def __init__(self, mu=0., s=1., *args, **kwargs):
+        super(Logistic, self).__init__(*args, **kwargs)
+
+        self.mu = tt.as_tensor_variable(mu)
+        self.s = tt.as_tensor_variable(s)
+
+        self.mean = self.mode = mu
+        self.variance = s**2 * np.pi**2 / 3.
+
+    def logp(self, value):
+        mu = self.mu
+        s = self.s
+
+        return bound(
+            -(value - mu) / s - tt.log(s) - 2 * tt.log1p(tt.exp(-(value - mu) / s)),
+            s > 0
+        )
+
+    def random(self, point=None, size=None, repeat=None):
+        mu, s = draw_values([self.mu, self.s], point=point)
+
+        return generate_samples(
+            stats.logistic.rvs,
+            loc=mu, scale=s,
+            dist_shape=self.shape,
+            size=size
+        )
+
+    def _repr_latex_(self, name=None, dist=None):
+        if dist is None:
+            dist = self
+        mu = dist.mu
+        s = dist.s
+        return r'${} \sim \text{{Logistic}}(\mathit{{mu}}={}, \mathit{{s}}={})$'.format(name,
+                                                                get_variable_name(mu),
+                                                                get_variable_name(s))
+
+
 class Interpolated(Continuous):
     R"""
     Univariate probability distribution defined as a linear interpolation

pymc3/tests/test_distributions.py

@@ -14,7 +14,7 @@
                              NegativeBinomial, Geometric, Exponential, ExGaussian, Normal,
                              Flat, LKJCorr, Wald, ChiSquared, HalfNormal, DiscreteUniform,
                              Bound, Uniform, Triangular, Binomial, SkewNormal, DiscreteWeibull,
-                             Gumbel, Interpolated, ZeroInflatedBinomial, HalfFlat, AR1)
+                             Gumbel, Logistic, Interpolated, ZeroInflatedBinomial, HalfFlat, AR1)
 from ..distributions import continuous
 from pymc3.theanof import floatX
 from numpy import array, inf, log, exp
@@ -895,6 +895,11 @@ def gumbel(value, mu, beta):
             return floatX(sp.gumbel_r.logpdf(value, loc=mu, scale=beta))
         self.pymc3_matches_scipy(Gumbel, R, {'mu': R, 'beta': Rplusbig}, gumbel)
 
+    def test_logistic(self):
+        self.pymc3_matches_scipy(Logistic, R, {'mu': R, 's': Rplus},
+                                 lambda value, mu, s: sp.logistic.logpdf(value, mu, s),
+                                 decimal=select_by_precision(float64=6, float32=1))
+
     def test_multidimensional_beta_construction(self):
         with Model():
             Beta('beta', alpha=1., beta=1., shape=(10, 20))