Logit-normal distribution

pymc3/distributions/continuous.py

@@ -17,6 +17,7 @@
 
 from .dist_math import bound, logpow, gammaln, betaln, std_cdf, i0, i1, alltrue_elemwise
 from .distribution import Continuous, draw_values, generate_samples, Bound
+from pymc3.math import logit
 
 __all__ = ['Uniform', 'Flat', 'Normal', 'Beta', 'Exponential', 'Laplace',
            'StudentT', 'Cauchy', 'HalfCauchy', 'Gamma', 'Weibull',
@@ -63,8 +64,7 @@ def assert_negative_support(var, label, distname, value=-1e-6):
 
 
 def get_tau_sd(tau=None, sd=None):
-    """
-    Find precision and standard deviation
+    """Find precision and standard deviation
 
     .. math::
         \tau = \frac{1}{\sigma^2}
@@ -104,8 +104,7 @@ def get_tau_sd(tau=None, sd=None):
 
 
 class Uniform(Continuous):
-    R"""
-    Continuous uniform log-likelihood.
+    R"""Continuous uniform log-likelihood.
 
     .. math::
 
@@ -153,8 +152,7 @@ def logp(self, value):
 
 
 class Flat(Continuous):
-    """
-    Uninformative log-likelihood that returns 0 regardless of
+    """Uninformative log-likelihood that returns 0 regardless of
     the passed value.
     """
 
@@ -170,8 +168,7 @@ def logp(self, value):
 
 
 class Normal(Continuous):
-    R"""
-    Univariate normal log-likelihood.
+    R"""Univariate normal log-likelihood.
 
     .. math::
 
@@ -233,8 +230,7 @@ def logp(self, value):
 
 
 class HalfNormal(PositiveContinuous):
-    R"""
-    Half-normal log-likelihood.
+    R"""Half-normal log-likelihood.
 
     .. math::
 
@@ -284,8 +280,7 @@ def logp(self, value):
 
 
 class Wald(PositiveContinuous):
-    R"""
-    Wald log-likelihood.
+    R"""Wald log-likelihood.
 
     .. math::
 
@@ -405,8 +400,7 @@ def logp(self, value):
 
 
 class Beta(UnitContinuous):
-    R"""
-    Beta log-likelihood.
+    R"""Beta log-likelihood.
 
     .. math::
 
@@ -493,8 +487,7 @@ def logp(self, value):
 
 
 class Exponential(PositiveContinuous):
-    R"""
-    Exponential log-likelihood.
+    R"""Exponential log-likelihood.
 
     .. math::
 
@@ -535,8 +528,7 @@ def logp(self, value):
 
 
 class Laplace(Continuous):
-    R"""
-    Laplace log-likelihood.
+    R"""Laplace log-likelihood.
 
     .. math::
 
@@ -580,8 +572,7 @@ def logp(self, value):
 
 
 class Lognormal(PositiveContinuous):
-    R"""
-    Log-normal log-likelihood.
+    R"""Log-normal log-likelihood.
 
     Distribution of any random variable whose logarithm is normally
     distributed. A variable might be modeled as log-normal if it can
@@ -597,7 +588,7 @@ class Lognormal(PositiveContinuous):
     ========  =========================================================================
     Support   :math:`x \in (0, 1)`
     Mean      :math:`\exp\{\mu + \frac{1}{2\tau}\}`
-    Variance  :math:\(\exp\{\frac{1}{\tau}\} - 1\) \times \exp\{2\mu + \frac{1}{\tau}\}
+    Variance  :math:`\(\exp\{\frac{1}{\tau}\} - 1\) \times \exp\{2\mu + \frac{1}{\tau}\}`
     ========  =========================================================================
 
     Parameters
@@ -643,9 +634,48 @@ def logp(self, value):
                      tau > 0)
 
 
+class Logitnormal(PositiveContinuous):
+    R"""Logit-normal log-likelihood.
+
+    .. math::
+
+       f(x \mid \mu, \tau) =
+           \frac{1}{x(x-1)} \sqrt{\frac{\tau}{2\pi}}
+           \exp\left\{ -\frac{\tau}{2} (\logit(x)-\mu)^2 \right\}
+
+    ========  =========================================================================
+    Support   :math:`x \in (0, 1)`
+    Mean      undefined
+    Variance  undefined
+    ========  =========================================================================
+
+    Parameters
+    ----------
+    mu : float
+        Location parameter.
+    tau : float
+        Scale parameter (tau > 0).
+    """
+
+    def __init__(self, mu=0.5, tau=None, *args, **kwargs):
+        super(Logitnormal, self).__init__(*args, **kwargs)
+
+        self.mu = mu = tt.as_tensor_variable(mu)
+        self.tau = tau = tt.as_tensor_variable(tau)
+
+        self.median = logit(mu)
+        assert_negative_support(tau, 'tau', 'Logitnormal')
+
+    def logp(self, value):
+        mu = self.mu
+        tau = self.tau
+        return bound(-0.5 * tau * (logit(value) - mu) ** 2
+                     + 0.5 * tt.log(tau / (2. * np.pi))
+                     - tt.log(value * (1 - value)), value > 0, value < 1, tau > 0)
+
+
 class StudentT(Continuous):
-    R"""
-    Non-central Student's T log-likelihood.
+    R"""Non-central Student's T log-likelihood.
 
     Describes a normal variable whose precision is gamma distributed.
     If only nu parameter is passed, this specifies a standard (central)
@@ -708,8 +738,7 @@ def logp(self, value):
 
 
 class Pareto(PositiveContinuous):
-    R"""
-    Pareto log-likelihood.
+    R"""Pareto log-likelihood.
 
     Often used to characterize wealth distribution, or other examples of the
     80/20 rule.
@@ -770,8 +799,7 @@ def logp(self, value):
 
 
 class Cauchy(Continuous):
-    R"""
-    Cauchy log-likelihood.
+    R"""Cauchy log-likelihood.
 
     Also known as the Lorentz or the Breit-Wigner distribution.
 
@@ -822,8 +850,7 @@ def logp(self, value):
 
 
 class HalfCauchy(PositiveContinuous):
-    R"""
-    Half-Cauchy log-likelihood.
+    R"""Half-Cauchy log-likelihood.
 
     .. math::
 
@@ -868,8 +895,7 @@ def logp(self, value):
 
 
 class Gamma(PositiveContinuous):
-    R"""
-    Gamma log-likelihood.
+    R"""Gamma log-likelihood.
 
     Represents the sum of alpha exponentially distributed random variables,
     each of which has mean beta.
@@ -951,8 +977,7 @@ def logp(self, value):
 
 
 class InverseGamma(PositiveContinuous):
-    R"""
-    Inverse gamma log-likelihood, the reciprocal of the gamma distribution.
+    R"""Inverse gamma log-likelihood, the reciprocal of the gamma distribution.
 
     .. math::
 
@@ -1012,8 +1037,7 @@ def logp(self, value):
 
 
 class ChiSquared(Gamma):
-    R"""
-    :math:`\chi^2` log-likelihood.
+    R""":math:`\chi^2` log-likelihood.
 
     .. math::
 
@@ -1038,8 +1062,7 @@ def __init__(self, nu, *args, **kwargs):
 
 
 class Weibull(PositiveContinuous):
-    R"""
-    Weibull log-likelihood.
+    R"""Weibull log-likelihood.
 
     .. math::
 
@@ -1102,8 +1125,7 @@ def StudentTpos(*args, **kwargs):
 
 
 class ExGaussian(Continuous):
-    R"""
-    Exponentially modified Gaussian log-likelihood.
+    R"""Exponentially modified Gaussian log-likelihood.
 
     Results from the convolution of a normal distribution with an exponential
     distribution.
@@ -1184,8 +1206,7 @@ def logp(self, value):
 
 
 class VonMises(Continuous):
-    R"""
-    Univariate VonMises log-likelihood.
+    R"""Univariate VonMises log-likelihood.
 
     .. math::
         f(x \mid \mu, \kappa) =
@@ -1232,8 +1253,7 @@ def logp(self, value):
 
 
 class SkewNormal(Continuous):
-    R"""
-    Univariate skew-normal log-likelihood.
+    R"""Univariate skew-normal log-likelihood.
 
     .. math::
        f(x \mid \mu, \tau, \alpha) =
@@ -1307,8 +1327,7 @@ def logp(self, value):
 
 
 class Triangular(Continuous):
-    """
-    Continuous Triangular log-likelihood
+    """Continuous Triangular log-likelihood
     Implemented by J. A. Fonseca 22/12/16
 
     Parameters