add skew-normal distribution and tests (#1392)

* add skew-normal distribution and tests

* replace the order of tau and sd to make it match the calling signature

Author: aloctavodia

pymc3/distributions/__init__.py

@@ -23,6 +23,7 @@
 from .continuous import InverseGamma
 from .continuous import ExGaussian
 from .continuous import VonMises
+from .continuous import SkewNormal
 
 from .discrete import Binomial
 from .discrete import BetaBinomial
@@ -110,5 +111,6 @@
            'LKJCorr',
            'AR1',
            'GaussianRandomWalk',
-           'GARCH11'
+           'GARCH11',
+           'SkewNormal'
            ]

pymc3/distributions/continuous.py

@@ -18,8 +18,9 @@
 
 __all__ = ['Uniform', 'Flat', 'Normal', 'Beta', 'Exponential', 'Laplace',
            'StudentT', 'Cauchy', 'HalfCauchy', 'Gamma', 'Weibull',
-           'Bound', 'HalfStudentT', 'StudentTpos', 'Lognormal', 'ChiSquared', 'HalfNormal',
-           'Wald', 'Pareto', 'InverseGamma', 'ExGaussian', 'VonMises']
+           'Bound', 'HalfStudentT', 'StudentTpos', 'Lognormal', 'ChiSquared',
+           'HalfNormal', 'Wald', 'Pareto', 'InverseGamma', 'ExGaussian',
+           'VonMises', 'SkewNormal']
 
 
 class PositiveContinuous(Continuous):
@@ -1237,3 +1238,66 @@ def logp(self, value):
         mu = self.mu
         kappa = self.kappa
         return bound(kappa * tt.cos(mu - value) - tt.log(2 * np.pi * i0(kappa)), value >= -np.pi, value <= np.pi, kappa >= 0)
+
+
+class SkewNormal(Continuous):
+    R"""
+    Univariate skew-normal log-likelihood.
+    .. math::
+       f(x \mid \mu, \tau, \alpha) = 
+       2 \Phi((x-\mu)\sqrt{\tau}\alpha) \phi(x,\mu,\tau)
+    ========  ==========================================
+    Support   :math:`x \in \mathbb{R}`
+    Mean      :math:`\mu + \sigma \sqrt{\frac{2}{\pi}} \frac {\alpha }{{\sqrt {1+\alpha ^{2}}}}`
+    Variance  :math:`\sigma^2 \left(  1-\frac{2\alpha^2}{(\alpha^2+1) \pi} \right)`
+    ========  ==========================================
+    Skew-normal distribution can be parameterized either in terms of precision
+    or standard deviation. The link between the two parametrizations is
+    given by
+    .. math::
+       \tau = \dfrac{1}{\sigma^2}
+    Parameters
+    ----------
+    mu : float
+        Location parameter.
+    sd : float
+        Scale parameter (sd > 0).
+    tau : float
+        Alternative scale parameter (tau > 0).
+    alpha : float
+        Skewness parameter.
+    
+    Notes
+    -----
+    When alpha=0 we recover the Normal distribution and mu becomes the mean, 
+    tau the precision and sd the standard deviation. In the limit of alpha 
+    approaching plus/minus infinite we get a half-normal distribution.       
+        
+    """
+    def __init__(self, mu=0.0, sd=None, tau=None, alpha=1,  *args, **kwargs):
+        super(SkewNormal, self).__init__(*args, **kwargs)
+        self.mu = mu
+        self.tau, self.sd = get_tau_sd(tau=tau, sd=sd)
+        self.alpha = alpha
+        self.mean = mu + self.sd * (2 / np.pi)**0.5 * alpha / (1 + alpha**2)**0.5
+        self.variance = self.sd**2 * (1 - (2 * alpha**2) / ((1 + alpha**2) * np.pi))
+
+    def random(self, point=None, size=None, repeat=None):
+        mu, tau, sd, alpha = draw_values(
+            [self.mu, self.tau, self.sd, self.alpha], point=point)
+        return generate_samples(stats.skewnorm.rvs,
+                                a=alpha, loc=mu, scale=tau**-0.5,
+                                dist_shape=self.shape,
+                                size=size)
+
+    def logp(self, value):
+        tau = self.tau
+        sd = self.sd
+        mu = self.mu
+        alpha = self.alpha
+        return bound(
+            tt.log(1 +
+            tt.erf(((value - mu) * tt.sqrt(tau) * alpha) / tt.sqrt(2)))
+            + (-tau * (value - mu)**2
+            + tt.log(tau / np.pi / 2.)) / 2.,
+            tau > 0, sd > 0)

pymc3/tests/test_distributions.py

@@ -12,7 +12,7 @@
                              Gamma, Cauchy, HalfCauchy, Lognormal, Laplace, NegativeBinomial,
                              Geometric, Exponential, ExGaussian, Normal, Flat, LKJCorr, Wald,
                              ChiSquared, HalfNormal, DiscreteUniform, Bound, Uniform,
-                             Binomial, Wishart)
+                             Binomial, Wishart, SkewNormal)
 from ..distributions import continuous
 from numpy import array, inf, log, exp
 from numpy.testing import assert_almost_equal
@@ -457,6 +457,10 @@ def test_student_tpos(self):
         self.pymc3_matches_scipy(HalfStudentT, Rplus, {'nu': Rplus, 'mu': R, 'lam': Rplus},
                                  lambda value, nu, mu, lam: sp.t.logpdf(value, nu, mu, lam**-.5))
 
+    def test_skew_normal(self):
+        self.pymc3_matches_scipy(SkewNormal, R, {'mu': R, 'sd': Rplusbig, 'alpha': R},
+                                 lambda value, alpha, mu, sd: sp.skewnorm.logpdf(value, alpha, mu, sd))
+
     def test_binomial(self):
         self.pymc3_matches_scipy(Binomial, Nat, {'n': NatSmall, 'p': Unit},
                                  lambda value, n, p: sp.binom.logpmf(value, n, p))

pymc3/tests/test_distributions_random.py

@@ -12,7 +12,7 @@
                              ZeroInflatedNegativeBinomial, ConstantDist, Poisson, Bernoulli, Beta,
                              BetaBinomial, StudentT, Weibull, Pareto, InverseGamma, Gamma, Cauchy,
                              HalfCauchy, Lognormal, Laplace, NegativeBinomial, Geometric,
-                             Exponential, ExGaussian, Normal, Flat, Wald, ChiSquared,
+                             Exponential, ExGaussian, Normal, SkewNormal, Flat, Wald, ChiSquared,
                              HalfNormal, DiscreteUniform, Bound, Uniform, Binomial, draw_values)
 from ..model import Model, Point
 
@@ -257,6 +257,9 @@ def check(self, dist, **kwargs):
 
     def test_normal(self):
         self.check(Normal, mu=0., tau=1.)
+        
+    def test_skew_normal(self):
+        self.check(SkewNormal, mu=0., sd=1., alpha=5.)
 
     def test_uniform(self):
         self.check(Uniform, lower=0., upper=1.)
@@ -359,6 +362,9 @@ def check(self, dist, **kwargs):
 
     def test_normal(self):
         self.check(Normal, mu=self.zeros, tau=self.ones)
+        
+    def test_skew_normal(self):
+        self.check(SkewNormal, mu=self.zeros, sd=self.ones, alpha=self.ones * 5)
 
     def test_uniform(self):
         self.check(Uniform, lower=self.zeros, upper=self.ones)
@@ -475,6 +481,9 @@ def check(self, dist, **kwargs):
 
     def test_normal(self):
         self.check(Normal, mu=self.zeros, tau=self.ones)
+        
+    def test_skew_normal(self):
+        self.check(SkewNormal, mu=self.zeros, sd=self.ones, alpha=self.ones * 5)
 
     def test_uniform(self):
         self.check(Uniform, lower=self.zeros, upper=self.ones)
@@ -593,6 +602,11 @@ def ref_rand(size, mu, sd):
             return st.norm.rvs(size=size, loc=mu, scale=sd)
         pymc3_random(Normal, {'mu': R, 'sd': Rplus}, ref_rand=ref_rand)
 
+    def test_skew_normal(self):
+        def ref_rand(size, alpha, mu, sd):
+            return st.skewnorm.rvs(size=size, a=alpha, loc=mu, scale=sd)
+        pymc3_random(SkewNormal, {'mu': R, 'sd': Rplus, 'alpha': R}, ref_rand=ref_rand)
+
     def test_half_normal(self):
         def ref_rand(size, tau):
             return st.halfnorm.rvs(size=size, loc=0, scale=tau ** -0.5)