ENH: added generalized normal and gen. half normal distributions to s…

…cipy.stats

THANKS.txt

@@ -34,7 +34,7 @@ Damian Eads for hierarchical clustering, dendrogram plotting,
 Anne Archibald for kd-trees and nearest neighbor in scipy.spatial.
 Pauli Virtanen for Sphinx documentation generation, online documentation
     framework and interpolation bugfixes.
-Josef Perktold for major improvements to scipy.stats and its test suite and 
+Josef Perktold for major improvements to scipy.stats and its test suite and
               fixes and tests to optimize.curve_fit and leastsq.
 David Morrill for getting the scoreboard test system up and running.
 Louis Luangkesorn for providing multiple tests for the stats module.
@@ -65,7 +65,7 @@ Ondrej Certik for Debian packaging.
 Paul Ivanov for porting Numeric-style C code to the new NumPy API.
 Ariel Rokem for contributions on percentileofscore fixes and tests.
 Yosef Meller for tests in the optimization module.
-Ralf Gommers for release management, code clean up and improvements 
+Ralf Gommers for release management, code clean up and improvements
     to doc-string generation.
 Bruce Southey for bug-fixes and improvements to scipy.stats.
 Ernest Adrogué for the Skellam distribution.
@@ -92,7 +92,7 @@ Thouis (Ray) Jones for bug fixes in ndimage.
 Yaroslav Halchenko for a bug fix in ndimage.
 Thomas Robitaille for the IDL 'save' reader.
 Fazlul Shahriar for fixes to the NetCDF3 I/O.
-Chris Jordan-Squire for bug fixes, documentation improvements and 
+Chris Jordan-Squire for bug fixes, documentation improvements and
     scipy.special.logit & expit.
 Christoph Gohlke for many bug fixes and help with Windows specific issues.
 Jacob Silterra for cwt-based peak finding in scipy.signal.
@@ -141,6 +141,7 @@ Florian Wilhelm for usage of RandomState in scipy.stats distributions.
 Robert T. McGibbon for Levinson-Durbin Toeplitz solver.
 Alex Conley for the Exponentially Modified Normal distribution.
 Abraham Escalante for contributions to scipy.stats
+Johannes Ballé for the generalized normal distribution.
 
 
 Institutions

scipy/stats/_continuous_distns.py

@@ -16,7 +16,7 @@
 
 from numpy import (where, arange, putmask, ravel, sum, shape,
                    log, sqrt, exp, arctanh, tan, sin, arcsin, arctan,
-                   tanh, cos, cosh, sinh)
+                   tanh, cos, cosh, sinh, sign)
 
 from numpy import polyval, place, extract, any, asarray, nan, inf, pi
 
@@ -4408,6 +4408,126 @@ def _entropy(self, c):
 wrapcauchy = wrapcauchy_gen(a=0.0, b=2*pi, name='wrapcauchy')
 
 
+class gennorm_gen(rv_continuous):
+    """A generalized normal continuous random variable.
+
+    %(before_notes)s
+
+    Notes
+    -----
+    The probability density function for `gennorm` is [1]_::
+
+                                     beta
+        gennorm.pdf(x, beta) =  ---------------  exp(-|x|**beta)
+                                2 gamma(1/beta)
+
+    `gennorm` takes ``beta`` as a shape parameter.
+    For ``beta = 1``, it is identical to a Laplace distribution.
+    For ``beta = 2``, it is identical to a normal distribution (with ``scale=1/sqrt(2)``).
+
+    See Also
+    --------
+    laplace: Laplace distribution
+    norm: normal distribution
+
+    References
+    ----------
+
+    .. [1] "Generalized normal distribution, Version 1",
+           https://en.wikipedia.org/wiki/Generalized_normal_distribution#Version_1
+
+    %(example)s
+
+    """
+
+    def _pdf(self, x, beta):
+        return exp(self._logpdf(x, beta))
+
+    def _logpdf(self, x, beta):
+        return log(.5 * beta) - special.gammaln(1. / beta) - abs(x)**beta
+
+    def _cdf(self, x, beta):
+        c = .5 * sign(x)
+        return (.5 + c) - c * special.gammaincc(1. / beta, abs(x)**beta)
+
+    def _ppf(self, x, beta):
+        c = sign(x - .5)
+        return c * special.gammainccinv(1. / beta, (1. + c) - 2. * c * x) ** (1. / beta)
+
+    def _sf(self, x, beta):
+        return self._cdf(-x, beta)
+
+    def _isf(self, x, beta):
+        return -self._ppf(x, beta)
+
+    def _stats(self, beta):
+        c1 = special.gammaln(1. / beta)
+        c3 = special.gammaln(3. / beta)
+        c5 = special.gammaln(5. / beta)
+        return 0., exp(c3 - c1), 0., exp(c5 + c1 - 2. * c3) - 3.
+
+    def _entropy(self, beta):
+        return 1. / beta - log(.5 * beta) + special.gammaln(1. / beta)
+
+gennorm = gennorm_gen(name='gennorm')
+
+class halfgennorm_gen(rv_continuous):
+    """The upper half of a generalized normal continuous random variable.
+
+    %(before_notes)s
+
+    Notes
+    -----
+    The probability density function for `halfgennorm` is::
+
+                                        beta
+        halfgennorm.pdf(x, beta) =  -------------  exp(-|x|**beta)
+                                    gamma(1/beta)
+
+    `gennorm` takes ``beta`` as a shape parameter.
+    For ``beta = 1``, it is identical to an exponential distribution.
+    For ``beta = 2``, it is identical to a half normal distribution (with ``scale=1/sqrt(2)``).
+
+    See Also
+    --------
+    gennorm: generalized normal distribution
+    expon: exponential distribution
+    halfnorm: half normal distribution
+
+    References
+    ----------
+
+    .. [1] "Generalized normal distribution, Version 1",
+           https://en.wikipedia.org/wiki/Generalized_normal_distribution#Version_1
+
+    %(example)s
+
+    """
+
+    def _pdf(self, x, beta):
+        return exp(self._logpdf(x, beta))
+
+    def _logpdf(self, x, beta):
+        return log(beta) - special.gammaln(1. / beta) - x**beta
+
+    def _cdf(self, x, beta):
+        return special.gammainc(1. / beta, x**beta)
+
+    def _ppf(self, x, beta):
+        return special.gammaincinv(1. / beta, x) ** (1. / beta)
+
+    def _sf(self, x, beta):
+        return special.gammaincc(1. / beta, x**beta)
+
+    def _isf(self, x, beta):
+        return special.gammainccinv(1. / beta, x) ** (1. / beta)
+
+    def _entropy(self, beta):
+        return 1. / beta - log(beta) + special.gammaln(1. / beta)
+
+halfgennorm = halfgennorm_gen(a=0, name='halfgennorm')
+
+
 # Collect names of classes and objects in this module.
 pairs = list(globals().items())
 _distn_names, _distn_gen_names = get_distribution_names(pairs, rv_continuous)

scipy/stats/_distr_params.py

@@ -36,6 +36,8 @@
     ['gengamma', (4.4162385429431925, 3.1193091679242761)],
     ['genhalflogistic', (0.77274727809929322,)],
     ['genlogistic', (0.41192440799679475,)],
+    ['gennorm', (1.2988442399460265,)],
+    ['halfgennorm', (0.6748054997000371,)],
     ['genpareto', (0.1,)],   # use case with finite moments
     ['gilbrat', ()],
     ['gompertz', (0.94743713075105251,)],

scipy/stats/tests/test_distributions.py

@@ -41,7 +41,7 @@
          'weibull_min','weibull_max','dweibull','maxwell','rayleigh',
          'genlogistic', 'logistic','gumbel_l','gumbel_r','gompertz',
          'hypsecant', 'laplace', 'reciprocal','triang','tukeylambda',
-         'vonmises', 'vonmises_line', 'pearson3']
+         'vonmises', 'vonmises_line', 'pearson3', 'gennorm', 'halfgennorm']
 
 
 def _assert_hasattr(a, b, msg=None):
@@ -288,6 +288,42 @@ def test_ppf(self):
         expected = array([1.0, 2.0, 3.0])
         assert_array_almost_equal(vals, expected)
 
+class TestGennorm(TestCase):
+    def test_laplace(self):
+        # test against Laplace (special case for beta=1)
+        points = [1, 2, 3]
+        pdf1 = stats.gennorm.pdf(points, 1)
+        pdf2 = stats.laplace.pdf(points)
+        assert_almost_equal(pdf1, pdf2)
+
+    def test_norm(self):
+        # test against normal (special case for beta=2)
+        points = [1, 2, 3]
+        pdf1 = stats.gennorm.pdf(points, 2)
+        pdf2 = stats.norm.pdf(points, scale=2**-.5)
+        assert_almost_equal(pdf1, pdf2)
+
+class TestHalfgennorm(TestCase):
+    def test_expon(self):
+        # test against exponential (special case for beta=1)
+        points = [1, 2, 3]
+        pdf1 = stats.halfgennorm.pdf(points, 1)
+        pdf2 = stats.expon.pdf(points)
+        assert_almost_equal(pdf1, pdf2)
+
+    def test_halfnorm(self):
+        # test against half normal (special case for beta=2)
+        points = [1, 2, 3]
+        pdf1 = stats.halfgennorm.pdf(points, 2)
+        pdf2 = stats.halfnorm.pdf(points, scale=2**-.5)
+        assert_almost_equal(pdf1, pdf2)
+
+    def test_gennorm(self):
+        # test against generalized normal
+        points = [1, 2, 3]
+        pdf1 = stats.halfgennorm.pdf(points, .497324)
+        pdf2 = stats.gennorm.pdf(points, .497324)
+        assert_almost_equal(pdf1, 2*pdf2)
 
 class TestTruncnorm(TestCase):
     def test_ppf_ticket1131(self):
@@ -775,7 +811,7 @@ def test_tail(self):  # Regression test for ticket 807
         assert_equal(stats.expon.cdf(1e-18), 1e-18)
         assert_equal(stats.expon.isf(stats.expon.sf(40)), 40)
 
-        
+
 class TestExponNorm(TestCase):
     def test_moments(self):
         # Some moment test cases based on non-loc/scaled formula
@@ -785,9 +821,9 @@ def get_moms(lam, sig, mu):
             opK2 = 1.0 + 1 / (lam*sig)**2
             exp_skew = 2 / (lam * sig)**3 * opK2**(-1.5)
             exp_kurt = 6.0 * (1 + (lam * sig)**2)**(-2)
-            return [mu + 1/lam, sig*sig + 1.0/(lam*lam), exp_skew, exp_kurt] 
+            return [mu + 1/lam, sig*sig + 1.0/(lam*lam), exp_skew, exp_kurt]
 
-        mu, sig, lam = 0, 1, 1        
+        mu, sig, lam = 0, 1, 1
         K = 1.0 / (lam * sig)
         sts = stats.exponnorm.stats(K, loc=mu, scale=sig, moments='mvsk')
         assert_almost_equal(sts, get_moms(lam, sig, mu))
@@ -803,7 +839,7 @@ def get_moms(lam, sig, mu):
         K = 1.0 / (lam * sig)
         sts = stats.exponnorm.stats(K, loc=mu, scale=sig, moments='mvsk')
         assert_almost_equal(sts, get_moms(lam, sig, mu))
-        
+
     def test_extremes_x(self):
         # Test for extreme values against overflows
         assert_almost_equal(stats.exponnorm.pdf(-900, 1), 0.0)
@@ -907,7 +943,7 @@ def test_logpdf(self):
         assert_allclose(b.pdf(x), np.exp(b.logpdf(x)))
 
     def test_cdf(self):
-        # regression test for gh-4030: Implementation of 
+        # regression test for gh-4030: Implementation of
         # scipy.stats.betaprime.cdf()
         x = stats.betaprime.cdf(0, 0.2, 0.3)
         assert_equal(x, 0.0)
@@ -916,7 +952,7 @@ def test_cdf(self):
         x = np.array([0.2, 0.5, 0.6])
         cdfs = stats.betaprime.cdf(x, alpha, beta)
         assert_(np.isfinite(cdfs).all())
-        
+
         # check the new cdf implementation vs generic one:
         gen_cdf = stats.rv_continuous._cdf_single
         cdfs_g = [gen_cdf(stats.betaprime, val, alpha, beta) for val in x]
@@ -2206,12 +2242,12 @@ def test_lomax_accuracy():
     # regression test for gh-4033
     p = stats.lomax.ppf(stats.lomax.cdf(1e-100,1),1)
     assert_allclose(p, 1e-100)
-    
+
 def test_gompertz_accuracy():
     # Regression test for gh-4031
     p = stats.gompertz.ppf(stats.gompertz.cdf(1e-100,1),1)
     assert_allclose(p, 1e-100)
-    
+
 def test_truncexpon_accuracy():
     # regression test for gh-4035
     p = stats.truncexpon.ppf(stats.truncexpon.cdf(1e-100,1),1)