[Scipy-svn] r6963 - in trunk/scipy/stats: . tests

scipy-svn@scip... scipy-svn@scip...
Sun Nov 28 07:37:29 CST 2010


Author: rgommers
Date: 2010-11-28 07:37:29 -0600 (Sun, 28 Nov 2010)
New Revision: 6963

Modified:
   trunk/scipy/stats/distributions.py
   trunk/scipy/stats/info.py
   trunk/scipy/stats/tests/test_continuous_basic.py
Log:
DEP: deprecate stats.invnorm, add a copy as invgauss. See #1172.

Modified: trunk/scipy/stats/distributions.py
===================================================================
--- trunk/scipy/stats/distributions.py	2010-11-28 13:37:08 UTC (rev 6962)
+++ trunk/scipy/stats/distributions.py	2010-11-28 13:37:29 UTC (rev 6963)
@@ -6,6 +6,7 @@
 #
 
 import math
+import warnings
 from copy import copy
 
 from scipy.misc import comb, derivative
@@ -57,10 +58,10 @@
            'weibull_max', 'genlogistic', 'genpareto', 'genexpon', 'genextreme',
            'gamma', 'gengamma', 'genhalflogistic', 'gompertz', 'gumbel_r',
            'gumbel_l', 'halfcauchy', 'halflogistic', 'halfnorm', 'hypsecant',
-           'gausshyper', 'invgamma', 'invnorm', 'invweibull', 'johnsonsb',
-           'johnsonsu', 'laplace', 'levy', 'levy_l', 'levy_stable',
-           'logistic', 'loggamma', 'loglaplace', 'lognorm', 'gilbrat',
-           'maxwell', 'mielke', 'nakagami', 'ncx2', 'ncf', 't',
+           'gausshyper', 'invgamma', 'invnorm', 'invgauss', 'invweibull',
+           'johnsonsb', 'johnsonsu', 'laplace', 'levy', 'levy_l',
+           'levy_stable', 'logistic', 'loggamma', 'loglaplace', 'lognorm',
+           'gilbrat', 'maxwell', 'mielke', 'nakagami', 'ncx2', 'ncf', 't',
            'nct', 'pareto', 'lomax', 'powerlaw', 'powerlognorm', 'powernorm',
            'rdist', 'rayleigh', 'reciprocal', 'rice', 'recipinvgauss',
            'semicircular', 'triang', 'truncexpon', 'truncnorm',
@@ -3247,30 +3248,67 @@
 ## Inverse Normal Distribution
 # scale is gamma from DATAPLOT and B from Regress
 
+_invnorm_msg = \
+"""The `invnorm` distribution will be renamed to `invgauss` after scipy 0.9"""
 class invnorm_gen(rv_continuous):
     def _rvs(self, mu):
+        warnings.warn(_invnorm_msg, DeprecationWarning)
         return mtrand.wald(mu, 1.0, size=self._size)
     def _pdf(self, x, mu):
+        warnings.warn(_invnorm_msg, DeprecationWarning)
         return 1.0/sqrt(2*pi*x**3.0)*exp(-1.0/(2*x)*((x-mu)/mu)**2)
     def _logpdf(self, x, mu):
+        warnings.warn(_invnorm_msg, DeprecationWarning)
         return -0.5*log(2*pi) - 1.5*log(x) - ((x-mu)/mu)**2/(2*x)
     def _cdf(self, x, mu):
+        warnings.warn(_invnorm_msg, DeprecationWarning)
         fac = sqrt(1.0/x)
         C1 = norm.cdf(fac*(x-mu)/mu)
         C1 += exp(2.0/mu)*norm.cdf(-fac*(x+mu)/mu)
         return C1
     def _stats(self, mu):
+        warnings.warn(_invnorm_msg, DeprecationWarning)
         return mu, mu**3.0, 3*sqrt(mu), 15*mu
 invnorm = invnorm_gen(a=0.0, name='invnorm', longname="An inverse normal",
                       shapes="mu",extradoc="""
 
 Inverse normal distribution
 
+NOTE: `invnorm` will be renamed to `invgauss` after scipy 0.9
+
 invnorm.pdf(x,mu) = 1/sqrt(2*pi*x**3) * exp(-(x-mu)**2/(2*x*mu**2))
 for x > 0.
 """
                       )
 
+## Inverse Gaussian Distribution (used to be called 'invnorm'
+# scale is gamma from DATAPLOT and B from Regress
+
+class invgauss_gen(rv_continuous):
+    def _rvs(self, mu):
+        return mtrand.wald(mu, 1.0, size=self._size)
+    def _pdf(self, x, mu):
+        return 1.0/sqrt(2*pi*x**3.0)*exp(-1.0/(2*x)*((x-mu)/mu)**2)
+    def _logpdf(self, x, mu):
+        return -0.5*log(2*pi) - 1.5*log(x) - ((x-mu)/mu)**2/(2*x)
+    def _cdf(self, x, mu):
+        fac = sqrt(1.0/x)
+        C1 = norm.cdf(fac*(x-mu)/mu)
+        C1 += exp(2.0/mu)*norm.cdf(-fac*(x+mu)/mu)
+        return C1
+    def _stats(self, mu):
+        return mu, mu**3.0, 3*sqrt(mu), 15*mu
+invgauss = invgauss_gen(a=0.0, name='invgauss', longname="An inverse Gaussian",
+                        shapes="mu",extradoc="""
+
+Inverse Gaussian distribution
+
+invgauss.pdf(x,mu) = 1/sqrt(2*pi*x**3) * exp(-(x-mu)**2/(2*x*mu**2))
+for x > 0.
+"""
+                      )
+
+
 ## Inverted Weibull
 
 class invweibull_gen(rv_continuous):
@@ -4136,7 +4174,7 @@
 
 # FIXME: PPF does not work.
 class recipinvgauss_gen(rv_continuous):
-    def _rvs(self, mu): #added, taken from invnorm
+    def _rvs(self, mu): #added, taken from invgauss
         return 1.0/mtrand.wald(mu, 1.0, size=self._size)
     def _pdf(self, x, mu):
         return 1.0/sqrt(2*pi*x)*exp(-(1-mu*x)**2.0 / (2*x*mu**2.0))
@@ -4403,7 +4441,7 @@
 
 ## Wald distribution (Inverse Normal with shape parameter mu=1.0)
 
-class wald_gen(invnorm_gen):
+class wald_gen(invgauss_gen):
     """A Wald continuous random variable.
 
     %(before_notes)s
@@ -4418,11 +4456,11 @@
     def _rvs(self):
         return mtrand.wald(1.0, 1.0, size=self._size)
     def _pdf(self, x):
-        return invnorm._pdf(x, 1.0)
+        return invgauss._pdf(x, 1.0)
     def _logpdf(self, x):
-        return invnorm._logpdf(x, 1.0)
+        return invgauss._logpdf(x, 1.0)
     def _cdf(self, x):
-        return invnorm._cdf(x, 1.0)
+        return invgauss._cdf(x, 1.0)
     def _stats(self):
         return 1.0, 1.0, 3.0, 15.0
 wald = wald_gen(a=0.0, name="wald", extradoc="""

Modified: trunk/scipy/stats/info.py
===================================================================
--- trunk/scipy/stats/info.py	2010-11-28 13:37:08 UTC (rev 6962)
+++ trunk/scipy/stats/info.py	2010-11-28 13:37:29 UTC (rev 6963)
@@ -81,6 +81,7 @@
    hypsecant         Hyperbolic Secant
    invgamma          Inverse Gamma
    invnorm           Inverse Normal
+   invgauss          Inverse Gaussian
    invweibull        Inverse Weibull
    johnsonsb         Johnson SB
    johnsonsu         Johnson SU

Modified: trunk/scipy/stats/tests/test_continuous_basic.py
===================================================================
--- trunk/scipy/stats/tests/test_continuous_basic.py	2010-11-28 13:37:08 UTC (rev 6962)
+++ trunk/scipy/stats/tests/test_continuous_basic.py	2010-11-28 13:37:29 UTC (rev 6963)
@@ -1,3 +1,5 @@
+import warnings
+
 import numpy.testing as npt
 import numpy as np
 import nose
@@ -15,7 +17,7 @@
 not for numerically exact results.
 
 
-TODO: 
+TODO:
 * make functioning test for skew and kurtosis
   still known failures - skip for now
 
@@ -70,6 +72,7 @@
     ['hypsecant', ()],
     ['invgamma', (2.0668996136993067,)],
     ['invnorm', (0.14546264555347513,)],
+    ['invgauss', (0.14546264555347513,)],
     ['invweibull', (10.58,)], # sample mean test fails at(0.58847112119264788,)]
     ['johnsonsb', (4.3172675099141058, 3.1837781130785063)],
     ['johnsonsu', (2.554395574161155, 2.2482281679651965)],
@@ -144,7 +147,7 @@
     'loglaplace', 'rdist', 'semicircular', 'invweibull', 'ksone',
     'cosine', 'kstwobign', 'truncnorm', 'mielke', 'recipinvgauss', 'levy',
     'johnsonsu', 'levy_l', 'powernorm', 'wrapcauchy',
-    'johnsonsb', 'truncexpon', 'rice', 'invnorm', 'invgamma',
+    'johnsonsb', 'truncexpon', 'rice', 'invnorm', 'invgauss', 'invgamma',
     'powerlognorm']
 
 distmiss = [[dist,args] for dist,args in distcont if dist in distmissing]
@@ -166,7 +169,7 @@
         skurt = stats.kurtosis(rvs)
         sskew = stats.skew(rvs)
         m,v = distfn.stats(*arg)
-        
+
         yield check_sample_meanvar_, distfn, arg, m, v, sm, sv, sn, distname + \
               'sample mean test'
         # the sample skew kurtosis test has known failures, not very good distance measure
@@ -219,14 +222,14 @@
     if not np.isinf(m):
         npt.assert_almost_equal(m1, m, decimal=10, err_msg= msg + \
                             ' - 1st moment')
-    else:                     # or np.isnan(m1), 
-        npt.assert_(np.isinf(m1), 
+    else:                     # or np.isnan(m1),
+        npt.assert_(np.isinf(m1),
                msg + ' - 1st moment -infinite, m1=%s' % str(m1))
         #np.isnan(m1) temporary special treatment for loggamma
     if not np.isinf(v):
         npt.assert_almost_equal(m2-m1*m1, v, decimal=10, err_msg= msg + \
                             ' - 2ndt moment')
-    else:                     #or np.isnan(m2), 
+    else:                     #or np.isnan(m2),
         npt.assert_(np.isinf(m2),
                msg + ' - 2nd moment -infinite, m2=%s' % str(m2))
         #np.isnan(m2) temporary special treatment for loggamma
@@ -356,6 +359,10 @@
         npt.assert_(pval > alpha, "D = " + str(D) + "; pval = " + str(pval) +
                "; alpha = " + str(alpha) + "\nargs = " + str(args))
 
+
+warnings.filterwarnings('ignore', message="The `invnorm` distribution")
+
+
 if __name__ == "__main__":
     #nose.run(argv=['', __file__])
     nose.runmodule(argv=[__file__,'-s'], exit=False)



More information about the Scipy-svn mailing list