diff --git a/pywafo/src/wafo/integrate.py b/pywafo/src/wafo/integrate.py index 0bd336a..35cbd10 100644 --- a/pywafo/src/wafo/integrate.py +++ b/pywafo/src/wafo/integrate.py @@ -272,9 +272,9 @@ def romberg(fun, a, b, releps=1e-3, abseps=1e-3): Ih2 = 0. Ih4 = rom[0, 0] abserr = Ih4 - #%epstab = zeros(1,decdigs+7) - #%newflg = 1 - #%[res,abserr,epstab,newflg] = dea(newflg,Ih4,abserr,epstab) + #epstab = zeros(1,decdigs+7) + #newflg = 1 + #[res,abserr,epstab,newflg] = dea(newflg,Ih4,abserr,epstab) two = 1 one = 0 for i in xrange(1, tableLimit): @@ -298,11 +298,11 @@ def romberg(fun, a, b, releps=1e-3, abseps=1e-3): if (2 <= i): [res, abserr] = dea3(Ih1, Ih2, Ih4) - #%Ih4 = res + # Ih4 = res if (abserr <= max(abseps, releps * abs(res))): break - #%rom(1,1:i) = rom(2,1:i) + # rom(1,1:i) = rom(2,1:i) two = one one = (one + 1) % 2 ipower *= 2 diff --git a/pywafo/src/wafo/source/rind2007/quad.mod b/pywafo/src/wafo/source/rind2007/quad.mod index 7068073..fa54374 100644 --- a/pywafo/src/wafo/source/rind2007/quad.mod +++ b/pywafo/src/wafo/source/rind2007/quad.mod @@ -1,5 +1,5 @@ -GFORTRAN module version '4' created from rind71mod.f on Sat May 05 23:15:45 2012 -MD5:0b1982321203177ab8efc0aefe21c275 -- If you edit this, you'll get what you deserve. +GFORTRAN module version '4' created from rind71mod.f on Fri Apr 05 14:43:36 2013 +MD5:7894124eb6e1c7a3f817b6596dd96b5e -- If you edit this, you'll get what you deserve. (() () () () () () () () () () () () () () () () () () () () () () () () () () ()) @@ -193,10 +193,11 @@ DIMENSION DUMMY) (REAL 8 0 0 REAL ()) 0 0 () (1 ASSUMED_SHAPE (CONSTANT DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0) 35 'xma' '' 'xma' 30 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0) -48 'xma' '' 'xma' 43 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 -DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0) -49 'n0' '' 'n0' 43 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY) +24 'n' '' 'n' 23 ((VARIABLE INOUT UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY) (INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0) +25 'wfout' '' 'wfout' 23 ((VARIABLE OUT UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 +DIMENSION DUMMY) (REAL 8 0 0 REAL ()) 0 0 () (1 ASSUMED_SHAPE (CONSTANT +(INTEGER 4 0 0 INTEGER ()) 0 '1') ()) 0 () () () 0 0) 26 'bpout' '' 'bpout' 23 ((VARIABLE OUT UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DIMENSION DUMMY) (REAL 8 0 0 REAL ()) 0 0 () (1 ASSUMED_SHAPE (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') ()) 0 () () () 0 0) @@ -220,11 +221,6 @@ DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0) DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0) 42 'n0' '' 'n0' 36 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY) (INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0) -24 'n' '' 'n' 23 ((VARIABLE INOUT UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY) -(INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0) -25 'wfout' '' 'wfout' 23 ((VARIABLE OUT UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 -DIMENSION DUMMY) (REAL 8 0 0 REAL ()) 0 0 () (1 ASSUMED_SHAPE (CONSTANT -(INTEGER 4 0 0 INTEGER ()) 0 '1') ()) 0 () () () 0 0) 10 'n' '' 'n' 9 ((VARIABLE INOUT UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY) (INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0) 11 'wfout' '' 'wfout' 9 ((VARIABLE OUT UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 @@ -253,16 +249,20 @@ DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0) DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0) 22 'n0' '' 'n0' 16 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY) (INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0) -44 'n' '' 'n' 43 ((VARIABLE INOUT UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY) -(INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0) 45 'wf' '' 'wf' 43 ((VARIABLE OUT UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DIMENSION DUMMY) (REAL 8 0 0 REAL ()) 0 0 () (1 ASSUMED_SHAPE (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') ()) 0 () () () 0 0) +44 'n' '' 'n' 43 ((VARIABLE INOUT UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY) +(INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0) 46 'bp' '' 'bp' 43 ((VARIABLE OUT UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DIMENSION DUMMY) (REAL 8 0 0 REAL ()) 0 0 () (1 ASSUMED_SHAPE (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') ()) 0 () () () 0 0) 47 'xmi' '' 'xmi' 43 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0) +48 'xma' '' 'xma' 43 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 +DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0) +49 'n0' '' 'n0' 43 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY) +(INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0) ) ('__convert_r4_r8' 0 8 'gausshe0' 0 2 'gaussla0' 0 7 'gaussle0' 0 3 diff --git a/pywafo/src/wafo/spectrum/core.py b/pywafo/src/wafo/spectrum/core.py index dfdf3ab..8ef0a1d 100644 --- a/pywafo/src/wafo/spectrum/core.py +++ b/pywafo/src/wafo/spectrum/core.py @@ -3259,7 +3259,7 @@ class SpecData2D(PlotData): pass def rotate(self, phi=0, rotateGrid=False, method='linear'): ''' - Rotate spectrum clockwise around the origin. + Rotate spectrum clockwise around the origin. Parameters ----------- diff --git a/pywafo/src/wafo/spectrum/test/test_models.py b/pywafo/src/wafo/spectrum/test/test_models.py index bf4a0d7..a3ad5d1 100644 --- a/pywafo/src/wafo/spectrum/test/test_models.py +++ b/pywafo/src/wafo/spectrum/test/test_models.py @@ -1,70 +1,76 @@ -import numpy as np -from wafo.spectrum.models import (Bretschneider, Jonswap, OchiHubble, Tmaspec, - Torsethaugen, McCormick, Wallop) - -def test_bretschneider(): - S = Bretschneider(Hm0=6.5,Tp=10) - vals = S((0,1,2,3)) - true_vals = np.array([ 0. , 1.69350993, 0.06352698, 0.00844783]) - assert((np.abs(vals-true_vals)<1e-7).all()) - -def test_if_jonswap_with_gamma_one_equals_bretschneider(): - S = Jonswap(Hm0=7, Tp=11,gamma=1) - vals = S((0,1,2,3)) - true_vals = np.array([ 0. , 1.42694133, 0.05051648, 0.00669692]) - assert((np.abs(vals-true_vals)<1e-7).all()) - - w = np.linspace(0,5) - S2 = Bretschneider(Hm0=7, Tp=11) - #JONSWAP with gamma=1 should be equal to Bretscneider: - assert(np.all(np.abs(S(w)-S2(w))<1.e-7)) - - -def test_tmaspec(): - S = Tmaspec(Hm0=7, Tp=11,gamma=1,h=10) - vals = S((0,1,2,3)) - true_vals = np.array([ 0. , 0.70106233, 0.05022433, 0.00669692]) - assert((np.abs(vals-true_vals)<1e-7).all()) - -def test_torsethaugen(): - - S = Torsethaugen(Hm0=7, Tp=11,gamma=1,h=10) - vals = S((0,1,2,3)) - true_vals = np.array([ 0. , 1.19989709, 0.05819794, 0.0093541 ]) - assert((np.abs(vals-true_vals)<1e-7).all()) - - vals = S.wind(range(4)) - true_vals = np.array([ 0. , 1.13560528, 0.05529849, 0.00888989]) - assert((np.abs(vals-true_vals)<1e-7).all()) - vals = S.swell(range(4)) - true_vals = np.array([ 0. , 0.0642918 , 0.00289946, 0.00046421]) - assert((np.abs(vals-true_vals)<1e-7).all()) - -def test_ochihubble(): - - S = OchiHubble(par=2) - vals = S(range(4)) - true_vals = np.array([ 0. , 0.90155636, 0.04185445, 0.00583207]) - assert((np.abs(vals-true_vals)<1e-7).all()) - -def test_mccormick(): - - S = McCormick(Hm0=6.5,Tp=10) - vals = S(range(4)) - true_vals = np.array([ 0. , 1.87865908, 0.15050447, 0.02994663]) - assert((np.abs(vals-true_vals)<1e-7).all()) - -def test_wallop(): - - S = Wallop(Hm0=6.5, Tp=10) - vals = S(range(4)) - true_vals = np.array([ 0.00000000e+00, 9.36921871e-01, 2.76991078e-03, - 7.72996150e-05]) - assert((np.abs(vals-true_vals)<1e-7).all()) - - -if __name__ == '__main__': - #main() - import nose -# nose.run() - test_tmaspec() \ No newline at end of file +import numpy as np +from wafo.spectrum.models import (Bretschneider, Jonswap, OchiHubble, Tmaspec, + Torsethaugen, McCormick, Wallop) + + +def test_bretschneider(): + S = Bretschneider(Hm0=6.5, Tp=10) + vals = S((0, 1, 2, 3)) + true_vals = np.array([0., 1.69350993, 0.06352698, 0.00844783]) + assert((np.abs(vals - true_vals) < 1e-7).all()) + + +def test_if_jonswap_with_gamma_one_equals_bretschneider(): + S = Jonswap(Hm0=7, Tp=11, gamma=1) + vals = S((0, 1, 2, 3)) + true_vals = np.array([0., 1.42694133, 0.05051648, 0.00669692]) + assert((np.abs(vals - true_vals) < 1e-7).all()) + + w = np.linspace(0, 5) + S2 = Bretschneider(Hm0=7, Tp=11) + # JONSWAP with gamma=1 should be equal to Bretscneider: + assert(np.all(np.abs(S(w) - S2(w)) < 1.e-7)) + + +def test_tmaspec(): + S = Tmaspec(Hm0=7, Tp=11, gamma=1, h=10) + vals = S((0, 1, 2, 3)) + true_vals = np.array([0., 0.70106233, 0.05022433, 0.00669692]) + assert((np.abs(vals - true_vals) < 1e-7).all()) + + +def test_torsethaugen(): + + S = Torsethaugen(Hm0=7, Tp=11, gamma=1, h=10) + vals = S((0, 1, 2, 3)) + true_vals = np.array([0., 1.19989709, 0.05819794, 0.0093541]) + assert((np.abs(vals - true_vals) < 1e-7).all()) + + vals = S.wind(range(4)) + true_vals = np.array([0., 1.13560528, 0.05529849, 0.00888989]) + assert((np.abs(vals - true_vals) < 1e-7).all()) + vals = S.swell(range(4)) + true_vals = np.array([0., 0.0642918, 0.00289946, 0.00046421]) + assert((np.abs(vals - true_vals) < 1e-7).all()) + + +def test_ochihubble(): + + S = OchiHubble(par=2) + vals = S(range(4)) + true_vals = np.array([0., 0.90155636, 0.04185445, 0.00583207]) + assert((np.abs(vals - true_vals) < 1e-7).all()) + + +def test_mccormick(): + + S = McCormick(Hm0=6.5, Tp=10) + vals = S(range(4)) + true_vals = np.array([0., 1.87865908, 0.15050447, 0.02994663]) + assert((np.abs(vals - true_vals) < 1e-7).all()) + + +def test_wallop(): + + S = Wallop(Hm0=6.5, Tp=10) + vals = S(range(4)) + true_vals = np.array([0.00000000e+00, 9.36921871e-01, 2.76991078e-03, + 7.72996150e-05]) + assert((np.abs(vals - true_vals) < 1e-7).all()) + + +if __name__ == '__main__': + # main() + import nose + nose.run() + #test_tmaspec() diff --git a/pywafo/src/wafo/spectrum/test/test_specdata1d.py b/pywafo/src/wafo/spectrum/test/test_specdata1d.py index 7626993..c9f7b17 100644 --- a/pywafo/src/wafo/spectrum/test/test_specdata1d.py +++ b/pywafo/src/wafo/spectrum/test/test_specdata1d.py @@ -1,203 +1,225 @@ -import wafo.spectrum.models as sm -from wafo.spectrum import SpecData1D -import numpy as np -def slow(f): - f.slow = True - return f - -@slow -def test_tocovmatrix(): - Sj = sm.Jonswap() - S = Sj.tospecdata() - acfmat = S.tocov_matrix(nr=3, nt=256, dt=0.1) - vals = acfmat[:2,:] - true_vals = np.array([[ 3.06073383, 0. , -1.67748256 , 0. ], - [ 3.05235423, -0.1674357 , -1.66811444, 0.18693242]]) - assert((np.abs(vals-true_vals)<1e-7).all()) - - -def test_tocovdata(): - Sj = sm.Jonswap() - S = Sj.tospecdata() - Nt = len(S.data)-1 - acf = S.tocovdata(nr=0, nt=Nt) - vals = acf.data[:5] - - true_vals = np.array([3.06090339, 2.22658399, 0.45307391, -1.17495501, -2.05649042]) - assert((np.abs(vals-true_vals)<1e-6).all()) - - -def test_to_t_pdf(): - ''' - The density of Tc is computed by: - >>> from wafo.spectrum import models as sm - >>> Sj = sm.Jonswap() - >>> S = Sj.tospecdata() - >>> f = S.to_t_pdf(pdef='Tc', paramt=(0, 10, 51), speed=7, seed=100) - >>> ['%2.3f' % val for val in f.data[:10]] - ['0.000', '0.014', '0.027', '0.040', '0.050', '0.059', '0.067', '0.072', '0.077', '0.081'] - - estimated error bounds - >>> ['%2.4f' % val for val in f.err[:10]] - ['0.0000', '0.0003', '0.0003', '0.0004', '0.0006', '0.0009', '0.0016', '0.0019', '0.0020', '0.0021'] - ''' -@slow -def test_sim(): - - Sj = sm.Jonswap();S = Sj.tospecdata() - ns =100; dt = .2 - x1 = S.sim(ns,dt=dt) - - - import scipy.stats as st - x2 = S.sim(20000,20) - truth1 = [0,np.sqrt(S.moment(1)[0]),0., 0.] - funs = [np.mean,np.std,st.skew,st.kurtosis] - for fun,trueval in zip(funs,truth1): - res = fun(x2[:,1::], axis=0) - m = res.mean() - sa = res.std() - #trueval, m, sa - assert(np.abs(m-trueval)>> import wafo.spectrum.models as sm - >>> import wafo.transform.models as wtm - >>> import wafo.objects as wo - >>> Hs = 7 - >>> Sj = sm.Jonswap(Hm0=Hs) - >>> S0 = Sj.tospecdata() - >>> ns =100; dt = .2 - >>> x1 = S0.sim(ns, dt=dt) - - >>> S = S0.copy() - >>> me, va, sk, ku = S.stats_nl(moments='mvsk') - >>> S.tr = wtm.TrHermite(mean=me, sigma=Hs/4, skew=sk, kurt=ku, ysigma=Hs/4) - >>> ys = wo.mat2timeseries(S.sim(ns=2**13)) - >>> g0, gemp = ys.trdata() - >>> t0 = g0.dist2gauss() - >>> t1 = S0.testgaussian(ns=2**13, t0=t0, cases=50) - >>> sum(t1>t0)<5 - True - ''' - -def test_moment(): - Sj = sm.Jonswap(Hm0=5) - S = Sj.tospecdata() #Make spectrum ob - vals, txt = S.moment() - true_vals = [1.5614600345079888, 0.95567089481941048] - true_txt = ['m0', 'm0tt'] - for tv,v in zip(true_vals, vals): - assert(tv==v) - -def test_nyquist_freq(): - - Sj = sm.Jonswap(Hm0=5) - S = Sj.tospecdata() #Make spectrum ob - assert(S.nyquist_freq()==3.0) - -def test_sampling_period(): - - Sj = sm.Jonswap(Hm0=5) - S = Sj.tospecdata() #Make spectrum ob - assert( S.sampling_period()== 1.0471975511965976) - -def test_normalize(): - - Sj = sm.Jonswap(Hm0=5) - S = Sj.tospecdata() #Make spectrum ob - S.moment(2) - ([1.5614600345079888, 0.95567089481941048], ['m0', 'm0tt']) - vals, txt = S.moment(2) - true_vals = [1.5614600345079888, 0.95567089481941048] - for tv,v in zip(true_vals, vals): - assert(tv==v) - - Sn = S.copy(); - Sn.normalize() - - #Now the moments should be one - new_vals, txt = Sn.moment(2) - for v in new_vals: - assert(np.abs(v-1.0)<1e-7) - -def test_characteristic(): - ''' - >>> import wafo.spectrum.models as sm - >>> Sj = sm.Jonswap(Hm0=5) - >>> S = Sj.tospecdata() #Make spectrum ob - >>> S.characteristic(1) - (array([ 8.59007646]), array([[ 0.03040216]]), ['Tm01']) - - >>> [ch, R, txt] = S.characteristic([1,2,3]) # fact a vector of integers - >>> ch; R; txt - array([ 8.59007646, 8.03139757, 5.62484314]) - array([[ 0.03040216, 0.02834263, nan], - [ 0.02834263, 0.0274645 , nan], - [ nan, nan, 0.01500249]]) - ['Tm01', 'Tm02', 'Tm24'] - - >>> S.characteristic('Ss') # fact a string - (array([ 0.04963112]), array([[ 2.63624782e-06]]), ['Ss']) - - >>> S.characteristic(['Hm0','Tm02']) # fact a list of strings - (array([ 4.99833578, 8.03139757]), array([[ 0.05292989, 0.02511371], - [ 0.02511371, 0.0274645 ]]), ['Hm0', 'Tm02']) - ''' -def test_bandwidth(): - - Sj = sm.Jonswap(Hm0=3, Tp=7) - w = np.linspace(0,4,256) - S = SpecData1D(Sj(w),w) #Make spectrum object from numerical values - vals = S.bandwidth([0,1,2,3]) - true_vals = np.array([ 0.73062845, 0.34476034, 0.68277527, 2.90817052]) - assert((np.abs(vals-true_vals)<1e-7).all()) - -def test_docstrings(): - import doctest - doctest.testmod() - -if __name__ == '__main__': - import nose - #nose.run() - #test_docstrings() - #test_tocovdata() - #test_tocovmatrix() - #test_sim() - #test_bandwidth() \ No newline at end of file +import wafo.spectrum.models as sm +from wafo.spectrum import SpecData1D +import numpy as np +import unittest + + +def slow(f): + f.slow = True + return f + + +class TestSpectrum(unittest.TestCase): + + @slow + def test_tocovmatrix(self): + Sj = sm.Jonswap() + S = Sj.tospecdata() + acfmat = S.tocov_matrix(nr=3, nt=256, dt=0.1) + vals = acfmat[:2, :] + true_vals = np.array([[3.06073383, 0.0000000, -1.67748256, 0.], + [3.05235423, -0.1674357, -1.66811444, 0.18693242]]) + self.assertTrue((np.abs(vals - true_vals) < 1e-7).all()) + + + +def test_tocovdata(): + Sj = sm.Jonswap() + S = Sj.tospecdata() + Nt = len(S.data) - 1 + acf = S.tocovdata(nr=0, nt=Nt) + vals = acf.data[:5] + + true_vals = np.array( + [3.06090339, 2.22658399, 0.45307391, -1.17495501, -2.05649042]) + assert((np.abs(vals - true_vals) < 1e-6).all()) + + +def test_to_t_pdf(): + Sj = sm.Jonswap() + S = Sj.tospecdata() + f = S.to_t_pdf(pdef='Tc', paramt=(0, 10, 51), speed=7, seed=100) + vals = ['%2.3f' % val for val in f.data[:10]] + truevals = ['0.000', '0.014', '0.027', '0.040', + '0.050', '0.059', '0.067', '0.072', '0.077', '0.081'] + + # estimated error bounds + vals = ['%2.4f' % val for val in f.err[:10]] + truevals = ['0.0000', '0.0003', '0.0003', '0.0004', + '0.0006', '0.0009', '0.0016', '0.0019', '0.0020', '0.0021'] + + +@slow +def test_sim(): + + Sj = sm.Jonswap() + S = Sj.tospecdata() + ns = 100 + dt = .2 + x1 = S.sim(ns, dt=dt) + + import scipy.stats as st + x2 = S.sim(20000, 20) + truth1 = [0, np.sqrt(S.moment(1)[0]), 0., 0.] + funs = [np.mean, np.std, st.skew, st.kurtosis] + for fun, trueval in zip(funs, truth1): + res = fun(x2[:, 1::], axis=0) + m = res.mean() + sa = res.std() + #trueval, m, sa + assert(np.abs(m - trueval) < sa) + + +@slow +def test_sim_nl(): + + Sj = sm.Jonswap() + S = Sj.tospecdata() + ns = 100 + dt = .2 + x1 = S.sim_nl(ns, dt=dt) + + import numpy as np + import scipy.stats as st + x2, x1 = S.sim_nl(ns=20000, cases=40) + truth1 = [0, np.sqrt(S.moment(1)[0][0])] + S.stats_nl(moments='sk') + truth1[-1] = truth1[-1] - 3 + + # truth1 + #[0, 1.7495200310090633, 0.18673120577479801, 0.061988521262417606] + + funs = [np.mean, np.std, st.skew, st.kurtosis] + for fun, trueval in zip(funs, truth1): + res = fun(x2[:, 1::], axis=0) + m = res.mean() + sa = res.std() + #trueval, m, sa + assert(np.abs(m - trueval) < 2 * sa) + + +def test_stats_nl(): + + Hs = 7. + Sj = sm.Jonswap(Hm0=Hs, Tp=11) + S = Sj.tospecdata() + me, va, sk, ku = S.stats_nl(moments='mvsk') + assert(me == 0.0) + assert(va == 3.0608203389019537) + assert(sk == 0.18673120577479801) + assert(ku == 3.0619885212624176) + + +def test_testgaussian(): + ''' + >>> import wafo.spectrum.models as sm + >>> import wafo.transform.models as wtm + >>> import wafo.objects as wo + >>> Hs = 7 + >>> Sj = sm.Jonswap(Hm0=Hs) + >>> S0 = Sj.tospecdata() + >>> ns =100; dt = .2 + >>> x1 = S0.sim(ns, dt=dt) + + >>> S = S0.copy() + >>> me, va, sk, ku = S.stats_nl(moments='mvsk') + >>> S.tr = wtm.TrHermite(mean=me, sigma=Hs/4, skew=sk, kurt=ku, ysigma=Hs/4) + >>> ys = wo.mat2timeseries(S.sim(ns=2**13)) + >>> g0, gemp = ys.trdata() + >>> t0 = g0.dist2gauss() + >>> t1 = S0.testgaussian(ns=2**13, t0=t0, cases=50) + >>> sum(t1>t0)<5 + True + ''' + + +def test_moment(): + Sj = sm.Jonswap(Hm0=5) + S = Sj.tospecdata() # Make spectrum ob + vals, txt = S.moment() + true_vals = [1.5614600345079888, 0.95567089481941048] + true_txt = ['m0', 'm0tt'] + for tv, v in zip(true_vals, vals): + assert(tv == v) + + +def test_nyquist_freq(): + + Sj = sm.Jonswap(Hm0=5) + S = Sj.tospecdata() # Make spectrum ob + assert(S.nyquist_freq() == 3.0) + + +def test_sampling_period(): + + Sj = sm.Jonswap(Hm0=5) + S = Sj.tospecdata() # Make spectrum ob + assert(S.sampling_period() == 1.0471975511965976) + + +def test_normalize(): + + Sj = sm.Jonswap(Hm0=5) + S = Sj.tospecdata() # Make spectrum ob + S.moment(2) + ([1.5614600345079888, 0.95567089481941048], ['m0', 'm0tt']) + vals, txt = S.moment(2) + true_vals = [1.5614600345079888, 0.95567089481941048] + for tv, v in zip(true_vals, vals): + assert(tv == v) + + Sn = S.copy() + Sn.normalize() + + # Now the moments should be one + new_vals, txt = Sn.moment(2) + for v in new_vals: + assert(np.abs(v - 1.0) < 1e-7) + + +def test_characteristic(): + ''' + >>> import wafo.spectrum.models as sm + >>> Sj = sm.Jonswap(Hm0=5) + >>> S = Sj.tospecdata() #Make spectrum ob + >>> S.characteristic(1) + (array([ 8.59007646]), array([[ 0.03040216]]), ['Tm01']) + + >>> [ch, R, txt] = S.characteristic([1,2,3]) # fact a vector of integers + >>> ch; R; txt + array([ 8.59007646, 8.03139757, 5.62484314]) + array([[ 0.03040216, 0.02834263, nan], + [ 0.02834263, 0.0274645 , nan], + [ nan, nan, 0.01500249]]) + ['Tm01', 'Tm02', 'Tm24'] + + >>> S.characteristic('Ss') # fact a string + (array([ 0.04963112]), array([[ 2.63624782e-06]]), ['Ss']) + + >>> S.characteristic(['Hm0','Tm02']) # fact a list of strings + (array([ 4.99833578, 8.03139757]), array([[ 0.05292989, 0.02511371], + [ 0.02511371, 0.0274645 ]]), ['Hm0', 'Tm02']) + ''' + + +def test_bandwidth(): + + Sj = sm.Jonswap(Hm0=3, Tp=7) + w = np.linspace(0, 4, 256) + S = SpecData1D(Sj(w), w) # Make spectrum object from numerical values + vals = S.bandwidth([0, 1, 2, 3]) + true_vals = np.array([0.73062845, 0.34476034, 0.68277527, 2.90817052]) + assert((np.abs(vals - true_vals) < 1e-7).all()) + + +def test_docstrings(): + import doctest + doctest.testmod() + +if __name__ == '__main__': + import nose + nose.run() + # test_docstrings() + # test_tocovdata() + # test_tocovmatrix() + # test_sim() + # test_bandwidth() diff --git a/pywafo/src/wafo/stats/six.py b/pywafo/src/wafo/stats/six.py index 09172d3..c537546 100644 --- a/pywafo/src/wafo/stats/six.py +++ b/pywafo/src/wafo/stats/six.py @@ -128,7 +128,6 @@ class MovedAttribute(_LazyDescr): return getattr(module, self.attr) - class _MovedItems(types.ModuleType): """Lazy loading of moved objects""" @@ -267,10 +266,12 @@ def iterkeys(d): """Return an iterator over the keys of a dictionary.""" return iter(getattr(d, _iterkeys)()) + def itervalues(d): """Return an iterator over the values of a dictionary.""" return iter(getattr(d, _itervalues)()) + def iteritems(d): """Return an iterator over the (key, value) pairs of a dictionary.""" return iter(getattr(d, _iteritems)()) @@ -279,8 +280,10 @@ def iteritems(d): if PY3: def b(s): return s.encode("latin-1") + def u(s): return s + if sys.version_info[1] <= 1: def int2byte(i): return bytes((i,)) @@ -293,6 +296,7 @@ if PY3: else: def b(s): return s + def u(s): return unicode(s, "unicode_escape") int2byte = chr @@ -306,13 +310,11 @@ if PY3: import builtins exec_ = getattr(builtins, "exec") - def reraise(tp, value, tb=None): if value.__traceback__ is not tb: raise value.with_traceback(tb) raise value - print_ = getattr(builtins, "print") del builtins @@ -329,17 +331,16 @@ else: locs = globs exec("""exec code in globs, locs""") - exec_("""def reraise(tp, value, tb=None): raise tp, value, tb """) - def print_(*args, **kwargs): """The new-style print function.""" fp = kwargs.pop("file", sys.stdout) if fp is None: return + def write(data): if not isinstance(data, basestring): data = str(data) @@ -385,4 +386,4 @@ _add_doc(reraise, """Reraise an exception.""") def with_metaclass(meta, base=object): """Create a base class with a metaclass.""" - return meta("NewBase", (base,), {}) \ No newline at end of file + return meta("NewBase", (base,), {}) diff --git a/pywafo/src/wafo/stats/tests/test_continuous_basic.py b/pywafo/src/wafo/stats/tests/test_continuous_basic.py new file mode 100644 index 0000000..7a9f71a --- /dev/null +++ b/pywafo/src/wafo/stats/tests/test_continuous_basic.py @@ -0,0 +1,390 @@ +from __future__ import division, print_function, absolute_import + +import warnings + +import numpy.testing as npt +import numpy as np +import nose + +from wafo import stats + +""" +Test all continuous distributions. + +Parameters were chosen for those distributions that pass the +Kolmogorov-Smirnov test. This provides safe parameters for each +distributions so that we can perform further testing of class methods. + +These tests currently check only/mostly for serious errors and exceptions, +not for numerically exact results. + + +TODO: +* make functioning test for skew and kurtosis + still known failures - skip for now + + +""" + +#currently not used +DECIMAL = 5 # specify the precision of the tests # increased from 0 to 5 +DECIMAL_kurt = 0 + +distcont = [ + ['alpha', (3.5704770516650459,)], + ['anglit', ()], + ['arcsine', ()], + ['beta', (2.3098496451481823, 0.62687954300963677)], + ['betaprime', (5, 6)], # avoid unbound error in entropy with (100, 86)], + ['bradford', (0.29891359763170633,)], + ['burr', (10.5, 4.3)], #incorrect mean and var for(0.94839838075366045, 4.3820284068855795)], + ['cauchy', ()], + ['chi', (78,)], + ['chi2', (55,)], + ['cosine', ()], + ['dgamma', (1.1023326088288166,)], + ['dweibull', (2.0685080649914673,)], + ['erlang', (20,)], #correction numargs = 1 + ['expon', ()], + ['exponpow', (2.697119160358469,)], + ['exponweib', (2.8923945291034436, 1.9505288745913174)], + ['f', (29, 18)], + ['fatiguelife', (29,)], #correction numargs = 1 + ['fisk', (3.0857548622253179,)], + ['foldcauchy', (4.7164673455831894,)], + ['foldnorm', (1.9521253373555869,)], + ['frechet_l', (3.6279911255583239,)], + ['frechet_r', (1.8928171603534227,)], + ['gamma', (1.9932305483800778,)], + ['gausshyper', (13.763771604130699, 3.1189636648681431, + 2.5145980350183019, 5.1811649903971615)], #veryslow + ['genexpon', (9.1325976465418908, 16.231956600590632, 3.2819552690843983)], + ['genextreme', (-0.1,)], # sample mean test fails for (3.3184017469423535,)], + ['gengamma', (4.4162385429431925, 3.1193091679242761)], + ['genhalflogistic', (0.77274727809929322,)], + ['genlogistic', (0.41192440799679475,)], + ['genpareto', (0.1,)], # use case with finite moments + ['gilbrat', ()], + ['gompertz', (0.94743713075105251,)], + ['gumbel_l', ()], + ['gumbel_r', ()], + ['halfcauchy', ()], + ['halflogistic', ()], + ['halfnorm', ()], + ['hypsecant', ()], + ['invgamma', (2.0668996136993067,)], + ['invgauss', (0.14546264555347513,)], + ['invweibull', (10.58,)], # sample mean test fails at(0.58847112119264788,)] + ['johnsonsb', (4.3172675099141058, 3.1837781130785063)], + ['johnsonsu', (2.554395574161155, 2.2482281679651965)], + ['ksone', (1000,)], #replace 22 by 100 to avoid failing range, ticket 956 + ['kstwobign', ()], + ['laplace', ()], + ['levy', ()], + ['levy_l', ()], +# ['levy_stable', (0.35667405469844993, +# -0.67450531578494011)], #NotImplementedError + # rvs not tested + ['loggamma', (0.41411931826052117,)], + ['logistic', ()], + ['loglaplace', (3.2505926592051435,)], + ['lognorm', (0.95368226960575331,)], + ['lomax', (1.8771398388773268,)], + ['maxwell', ()], + ['mielke', (10.4, 3.6)], # sample mean test fails for (4.6420495492121487, 0.59707419545516938)], + # mielke: good results if 2nd parameter >2, weird mean or var below + ['nakagami', (4.9673794866666237,)], + ['ncf', (27, 27, 0.41578441799226107)], + ['nct', (14, 0.24045031331198066)], + ['ncx2', (21, 1.0560465975116415)], + ['norm', ()], + ['pareto', (2.621716532144454,)], + ['pearson3', (0.1,)], + ['powerlaw', (1.6591133289905851,)], + ['powerlognorm', (2.1413923530064087, 0.44639540782048337)], + ['powernorm', (4.4453652254590779,)], + ['rayleigh', ()], + ['rdist', (0.9,)], # feels also slow +# ['rdist', (3.8266985793976525,)], #veryslow, especially rvs + #['rdist', (541.0,)], # from ticket #758 #veryslow + ['recipinvgauss', (0.63004267809369119,)], + ['reciprocal', (0.0062309367010521255, 1.0062309367010522)], + ['rice', (0.7749725210111873,)], + ['semicircular', ()], + ['t', (2.7433514990818093,)], + ['triang', (0.15785029824528218,)], + ['truncexpon', (4.6907725456810478,)], + ['truncnorm', (-1.0978730080013919, 2.7306754109031979)], + ['tukeylambda', (3.1321477856738267,)], + ['uniform', ()], + ['vonmises', (3.9939042581071398,)], + ['wald', ()], + ['weibull_max', (2.8687961709100187,)], + ['weibull_min', (1.7866166930421596,)], + ['wrapcauchy', (0.031071279018614728,)]] + +# for testing only specific functions +##distcont = [ +## ['erlang', (20,)], #correction numargs = 1 +## ['fatiguelife', (29,)], #correction numargs = 1 +## ['loggamma', (0.41411931826052117,)]] + +# for testing ticket:767 +##distcont = [ +## ['genextreme', (3.3184017469423535,)], +## ['genextreme', (0.01,)], +## ['genextreme', (0.00001,)], +## ['genextreme', (0.0,)], +## ['genextreme', (-0.01,)] +## ] + +##distcont = [['gumbel_l', ()], +## ['gumbel_r', ()], +## ['norm', ()] +## ] + +##distcont = [['norm', ()]] + +distmissing = ['wald', 'gausshyper', 'genexpon', 'rv_continuous', + 'loglaplace', 'rdist', 'semicircular', 'invweibull', 'ksone', + 'cosine', 'kstwobign', 'truncnorm', 'mielke', 'recipinvgauss', 'levy', + 'johnsonsu', 'levy_l', 'powernorm', 'wrapcauchy', + 'johnsonsb', 'truncexpon', 'rice', 'invgauss', 'invgamma', + 'powerlognorm'] + +distmiss = [[dist,args] for dist,args in distcont if dist in distmissing] +distslow = ['rdist', 'gausshyper', 'recipinvgauss', 'ksone', 'genexpon', + 'vonmises', 'rice', 'mielke', 'semicircular', 'cosine', 'invweibull', + 'powerlognorm', 'johnsonsu', 'kstwobign'] +#distslow are sorted by speed (very slow to slow) + +def _silence_fp_errors(func): + def wrap(*a, **kw): + olderr = np.seterr(all='ignore') + try: + return func(*a, **kw) + finally: + np.seterr(**olderr) + wrap.__name__ = func.__name__ + return wrap + +@_silence_fp_errors +def test_cont_basic(): + # this test skips slow distributions + for distname, arg in distcont[:]: + if distname in distslow: + continue + distfn = getattr(stats, distname) + np.random.seed(765456) + sn = 1000 + rvs = distfn.rvs(size=sn,*arg) + sm = rvs.mean() + sv = rvs.var() + 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 + #yield check_sample_skew_kurt, distfn, arg, sskew, skurt, distname + yield check_moment, distfn, arg, m, v, distname + yield check_cdf_ppf, distfn, arg, distname + yield check_sf_isf, distfn, arg, distname + yield check_pdf, distfn, arg, distname + if distname in ['wald']: + continue + yield check_pdf_logpdf, distfn, arg, distname + yield check_cdf_logcdf, distfn, arg, distname + yield check_sf_logsf, distfn, arg, distname + if distname in distmissing: + alpha = 0.01 + yield check_distribution_rvs, distname, arg, alpha, rvs + +@npt.dec.slow +def test_cont_basic_slow(): + # same as above for slow distributions + for distname, arg in distcont[:]: + if distname not in distslow: continue + distfn = getattr(stats, distname) + np.random.seed(765456) + sn = 1000 + rvs = distfn.rvs(size=sn,*arg) + sm = rvs.mean() + sv = rvs.var() + 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 + #yield check_sample_skew_kurt, distfn, arg, sskew, skurt, distname + yield check_moment, distfn, arg, m, v, distname + yield check_cdf_ppf, distfn, arg, distname + yield check_sf_isf, distfn, arg, distname + yield check_pdf, distfn, arg, distname + yield check_pdf_logpdf, distfn, arg, distname + yield check_cdf_logcdf, distfn, arg, distname + yield check_sf_logsf, distfn, arg, distname + #yield check_oth, distfn, arg # is still missing + if distname in distmissing: + alpha = 0.01 + yield check_distribution_rvs, distname, arg, alpha, rvs + +@_silence_fp_errors +def check_moment(distfn, arg, m, v, msg): + m1 = distfn.moment(1,*arg) + m2 = distfn.moment(2,*arg) + 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), + 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), + npt.assert_(np.isinf(m2), + msg + ' - 2nd moment -infinite, m2=%s' % str(m2)) + #np.isnan(m2) temporary special treatment for loggamma + +@_silence_fp_errors +def check_sample_meanvar_(distfn, arg, m, v, sm, sv, sn, msg): + #this did not work, skipped silently by nose + #check_sample_meanvar, sm, m, msg + 'sample mean test' + #check_sample_meanvar, sv, v, msg + 'sample var test' + if not np.isinf(m): + check_sample_mean(sm, sv, sn, m) + if not np.isinf(v): + check_sample_var(sv, sn, v) +## check_sample_meanvar( sm, m, msg + 'sample mean test') +## check_sample_meanvar( sv, v, msg + 'sample var test') + +def check_sample_mean(sm,v,n, popmean): + """ +from stats.stats.ttest_1samp(a, popmean): +Calculates the t-obtained for the independent samples T-test on ONE group +of scores a, given a population mean. + +Returns: t-value, two-tailed prob +""" +## a = asarray(a) +## x = np.mean(a) +## v = np.var(a, ddof=1) +## n = len(a) + df = n-1 + svar = ((n-1)*v) / float(df) #looks redundant + t = (sm-popmean)/np.sqrt(svar*(1.0/n)) + prob = stats.betai(0.5*df,0.5,df/(df+t*t)) + + #return t,prob + npt.assert_(prob > 0.01, 'mean fail, t,prob = %f, %f, m,sm=%f,%f' % (t,prob,popmean,sm)) + +def check_sample_var(sv,n, popvar): + ''' +two-sided chisquare test for sample variance equal to hypothesized variance + ''' + df = n-1 + chi2 = (n-1)*popvar/float(popvar) + pval = stats.chisqprob(chi2,df)*2 + npt.assert_(pval > 0.01, 'var fail, t,pval = %f, %f, v,sv=%f,%f' % (chi2,pval,popvar,sv)) + + +def check_sample_skew_kurt(distfn, arg, ss, sk, msg): + skew,kurt = distfn.stats(moments='sk',*arg) +## skew = distfn.stats(moment='s',*arg)[()] +## kurt = distfn.stats(moment='k',*arg)[()] + check_sample_meanvar( sk, kurt, msg + 'sample kurtosis test') + check_sample_meanvar( ss, skew, msg + 'sample skew test') + +def check_sample_meanvar(sm,m,msg): + if not np.isinf(m) and not np.isnan(m): + npt.assert_almost_equal(sm, m, decimal=DECIMAL, err_msg= msg + \ + ' - finite moment') +## else: +## npt.assert_(abs(sm) > 10000), msg='infinite moment, sm = ' + str(sm)) + +@_silence_fp_errors +def check_cdf_ppf(distfn,arg,msg): + values = [0.001, 0.5, 0.999] + npt.assert_almost_equal(distfn.cdf(distfn.ppf(values, *arg), *arg), + values, decimal=DECIMAL, err_msg= msg + \ + ' - cdf-ppf roundtrip') + +@_silence_fp_errors +def check_sf_isf(distfn,arg,msg): + npt.assert_almost_equal(distfn.sf(distfn.isf([0.1,0.5,0.9], *arg), *arg), + [0.1,0.5,0.9], decimal=DECIMAL, err_msg= msg + \ + ' - sf-isf roundtrip') + npt.assert_almost_equal(distfn.cdf([0.1,0.9], *arg), + 1.0-distfn.sf([0.1,0.9], *arg), + decimal=DECIMAL, err_msg= msg + \ + ' - cdf-sf relationship') + +@_silence_fp_errors +def check_pdf(distfn, arg, msg): + # compares pdf at median with numerical derivative of cdf + median = distfn.ppf(0.5, *arg) + eps = 1e-6 + pdfv = distfn.pdf(median, *arg) + if (pdfv < 1e-4) or (pdfv > 1e4): + # avoid checking a case where pdf is close to zero or huge (singularity) + median = median + 0.1 + pdfv = distfn.pdf(median, *arg) + cdfdiff = (distfn.cdf(median + eps, *arg) - + distfn.cdf(median - eps, *arg))/eps/2.0 + #replace with better diff and better test (more points), + #actually, this works pretty well + npt.assert_almost_equal(pdfv, cdfdiff, + decimal=DECIMAL, err_msg= msg + ' - cdf-pdf relationship') + +@_silence_fp_errors +def check_pdf_logpdf(distfn, args, msg): + # compares pdf at several points with the log of the pdf + points = np.array([0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8]) + vals = distfn.ppf(points, *args) + pdf = distfn.pdf(vals, *args) + logpdf = distfn.logpdf(vals, *args) + pdf = pdf[pdf != 0] + logpdf = logpdf[np.isfinite(logpdf)] + npt.assert_almost_equal(np.log(pdf), logpdf, decimal=7, err_msg=msg + " - logpdf-log(pdf) relationship") + +@_silence_fp_errors +def check_sf_logsf(distfn, args, msg): + # compares sf at several points with the log of the sf + points = np.array([0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8]) + vals = distfn.ppf(points, *args) + sf = distfn.sf(vals, *args) + logsf = distfn.logsf(vals, *args) + sf = sf[sf != 0] + logsf = logsf[np.isfinite(logsf)] + npt.assert_almost_equal(np.log(sf), logsf, decimal=7, err_msg=msg + " - logsf-log(sf) relationship") + +@_silence_fp_errors +def check_cdf_logcdf(distfn, args, msg): + # compares cdf at several points with the log of the cdf + points = np.array([0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8]) + vals = distfn.ppf(points, *args) + cdf = distfn.cdf(vals, *args) + logcdf = distfn.logcdf(vals, *args) + cdf = cdf[cdf != 0] + logcdf = logcdf[np.isfinite(logcdf)] + npt.assert_almost_equal(np.log(cdf), logcdf, decimal=7, err_msg=msg + " - logcdf-log(cdf) relationship") + + +@_silence_fp_errors +def check_distribution_rvs(dist, args, alpha, rvs): + #test from scipy.stats.tests + #this version reuses existing random variables + D,pval = stats.kstest(rvs, dist, args=args, N=1000) + if (pval < alpha): + D,pval = stats.kstest(dist,'',args=args, N=1000) + npt.assert_(pval > alpha, "D = " + str(D) + "; pval = " + str(pval) + + "; alpha = " + str(alpha) + "\nargs = " + str(args)) + + +if __name__ == "__main__": + #nose.run(argv=['', __file__]) + nose.runmodule(argv=[__file__,'-s'], exit=False) diff --git a/pywafo/src/wafo/stats/tests/test_continuous_extra.py b/pywafo/src/wafo/stats/tests/test_continuous_extra.py new file mode 100644 index 0000000..77c662f --- /dev/null +++ b/pywafo/src/wafo/stats/tests/test_continuous_extra.py @@ -0,0 +1,102 @@ +# contains additional tests for continuous distributions +# +# NOTE: one test, _est_cont_skip, that is renamed so that nose doesn't +# run it, +# 6 distributions return nan for entropy +# truncnorm fails by design for private method _ppf test +from __future__ import division, print_function, absolute_import + + +import numpy.testing as npt +import numpy as np + +from wafo import stats + +from test_continuous_basic import distcont + +DECIMAL = 5 + +@npt.dec.slow +def test_cont_extra(): + for distname, arg in distcont[:]: + distfn = getattr(stats, distname) + + yield check_ppf_limits, distfn, arg, distname + \ + ' ppf limit test' + yield check_isf_limits, distfn, arg, distname + \ + ' isf limit test' + yield check_loc_scale, distfn, arg, distname + \ + ' loc, scale test' + +@npt.dec.slow +def _est_cont_skip(): + for distname, arg in distcont: + distfn = getattr(stats, distname) + #entropy test checks only for isnan, currently 6 isnan left + yield check_entropy, distfn, arg, distname + \ + ' entropy nan test' + # _ppf test has 1 failure be design + yield check_ppf_private, distfn, arg, distname + \ + ' _ppf private test' + +def test_540_567(): + # test for nan returned in tickets 540, 567 + npt.assert_almost_equal(stats.norm.cdf(-1.7624320982),0.03899815971089126, + decimal=10, err_msg = 'test_540_567') + npt.assert_almost_equal(stats.norm.cdf(-1.7624320983),0.038998159702449846, + decimal=10, err_msg = 'test_540_567') + npt.assert_almost_equal(stats.norm.cdf(1.38629436112, loc=0.950273420309, + scale=0.204423758009),0.98353464004309321, + decimal=10, err_msg = 'test_540_567') + + +def check_ppf_limits(distfn,arg,msg): + below,low,upp,above = distfn.ppf([-1,0,1,2], *arg) + #print distfn.name, distfn.a, low, distfn.b, upp + #print distfn.name,below,low,upp,above + assert_equal_inf_nan(distfn.a,low, msg + 'ppf lower bound') + assert_equal_inf_nan(distfn.b,upp, msg + 'ppf upper bound') + npt.assert_(np.isnan(below), msg + 'ppf out of bounds - below') + npt.assert_(np.isnan(above), msg + 'ppf out of bounds - above') + +def check_ppf_private(distfn,arg,msg): + #fails by design for trunk norm self.nb not defined + ppfs = distfn._ppf(np.array([0.1,0.5,0.9]), *arg) + npt.assert_(not np.any(np.isnan(ppfs)), msg + 'ppf private is nan') + + +def check_isf_limits(distfn,arg,msg): + below,low,upp,above = distfn.isf([-1,0,1,2], *arg) + #print distfn.name, distfn.a, low, distfn.b, upp + #print distfn.name,below,low,upp,above + assert_equal_inf_nan(distfn.a,upp, msg + 'isf lower bound') + assert_equal_inf_nan(distfn.b,low, msg + 'isf upper bound') + npt.assert_(np.isnan(below), msg + 'isf out of bounds - below') + npt.assert_(np.isnan(above), msg + 'isf out of bounds - above') + + +def check_loc_scale(distfn,arg,msg): + m,v = distfn.stats(*arg) + loc, scale = 10.0, 10.0 + mt,vt = distfn.stats(loc=loc, scale=scale, *arg) + assert_equal_inf_nan(m*scale+loc,mt,msg + 'mean') + assert_equal_inf_nan(v*scale*scale,vt,msg + 'var') + +def check_entropy(distfn,arg,msg): + ent = distfn.entropy(*arg) + #print 'Entropy =', ent + npt.assert_(not np.isnan(ent), msg + 'test Entropy is nan') + +def assert_equal_inf_nan(v1,v2,msg): + npt.assert_(not np.isnan(v1)) + if not np.isinf(v1): + npt.assert_almost_equal(v1, v2, decimal=DECIMAL, err_msg = msg + \ + ' - finite') + else: + npt.assert_(np.isinf(v2) or np.isnan(v2), + msg + ' - infinite, v2=%s' % str(v2)) + +if __name__ == "__main__": + import nose + #nose.run(argv=['', __file__]) + nose.runmodule(argv=[__file__,'-s'], exit=False) diff --git a/pywafo/src/wafo/stats/tests/test_discrete_basic.py b/pywafo/src/wafo/stats/tests/test_discrete_basic.py new file mode 100644 index 0000000..659f809 --- /dev/null +++ b/pywafo/src/wafo/stats/tests/test_discrete_basic.py @@ -0,0 +1,268 @@ +import numpy.testing as npt +import numpy as np +import nose + +from wafo import stats + +DECIMAL_meanvar = 0#1 # was 0 + +distdiscrete = [ + ['bernoulli',(0.3,)], + ['binom', (5, 0.4)], + ['boltzmann',(1.4, 19)], + ['dlaplace', (0.8,)], #0.5 + ['geom', (0.5,)], + ['hypergeom',(30, 12, 6)], + ['hypergeom',(21,3,12)], #numpy.random (3,18,12) numpy ticket:921 + ['hypergeom',(21,18,11)], #numpy.random (18,3,11) numpy ticket:921 + ['logser', (0.6,)], # reenabled, numpy ticket:921 + ['nbinom', (5, 0.5)], + ['nbinom', (0.4, 0.4)], #from tickets: 583 + ['planck', (0.51,)], #4.1 + ['poisson', (0.6,)], + ['randint', (7, 31)], + ['skellam', (15, 8)]] +# ['zipf', (4,)] ] # arg=4 is ok, + # Zipf broken for arg = 2, e.g. weird .stats + # looking closer, mean, var should be inf for arg=2 + + +#@npt.dec.slow +def test_discrete_basic(): + for distname, arg in distdiscrete: + distfn = getattr(stats,distname) + #assert stats.dlaplace.rvs(0.8) is not None + np.random.seed(9765456) + rvs = distfn.rvs(size=2000,*arg) + supp = np.unique(rvs) + m,v = distfn.stats(*arg) + #yield npt.assert_almost_equal(rvs.mean(), m, decimal=4,err_msg='mean') + #yield npt.assert_almost_equal, rvs.mean(), m, 2, 'mean' # does not work + yield check_sample_meanvar, rvs.mean(), m, distname + ' sample mean test' + yield check_sample_meanvar, rvs.var(), v, distname + ' sample var test' + yield check_cdf_ppf, distfn, arg, distname + ' cdf_ppf' + yield check_cdf_ppf2, distfn, arg, supp, distname + ' cdf_ppf' + yield check_pmf_cdf, distfn, arg, distname + ' pmf_cdf' + + # zipf doesn't fail, but generates floating point warnings. + # Should be checked. + if not distname in ['zipf']: + yield check_oth, distfn, arg, distname + ' oth' + skurt = stats.kurtosis(rvs) + sskew = stats.skew(rvs) + yield check_sample_skew_kurt, distfn, arg, skurt, sskew, \ + distname + ' skew_kurt' + + # dlaplace doesn't fail, but generates lots of floating point warnings. + # Should be checked. + if not distname in ['dlaplace']: #['logser']: #known failure, fixed + alpha = 0.01 + yield check_discrete_chisquare, distfn, arg, rvs, alpha, \ + distname + ' chisquare' + +@npt.dec.slow +def test_discrete_extra(): + for distname, arg in distdiscrete: + distfn = getattr(stats,distname) + yield check_ppf_limits, distfn, arg, distname + \ + ' ppf limit test' + yield check_isf_limits, distfn, arg, distname + \ + ' isf limit test' + yield check_entropy, distfn, arg, distname + \ + ' entropy nan test' + +@npt.dec.skipif(True) +def test_discrete_private(): + #testing private methods mostly for debugging + # some tests might fail by design, + # e.g. incorrect definition of distfn.a and distfn.b + for distname, arg in distdiscrete: + distfn = getattr(stats,distname) + rvs = distfn.rvs(size=10000,*arg) + m,v = distfn.stats(*arg) + + yield check_ppf_ppf, distfn, arg + yield check_cdf_ppf_private, distfn, arg, distname + yield check_generic_moment, distfn, arg, m, 1, 3 # last is decimal + yield check_generic_moment, distfn, arg, v+m*m, 2, 3 # last is decimal + yield check_moment_frozen, distfn, arg, m, 1, 3 # last is decimal + yield check_moment_frozen, distfn, arg, v+m*m, 2, 3 # last is decimal + + +def check_sample_meanvar(sm,m,msg): + if not np.isinf(m): + npt.assert_almost_equal(sm, m, decimal=DECIMAL_meanvar, err_msg=msg + \ + ' - finite moment') + else: + npt.assert_(sm > 10000, msg='infinite moment, sm = ' + str(sm)) + +def check_sample_var(sm,m,msg): + npt.assert_almost_equal(sm, m, decimal=DECIMAL_meanvar, err_msg= msg + 'var') + +def check_cdf_ppf(distfn,arg,msg): + ppf05 = distfn.ppf(0.5,*arg) + cdf05 = distfn.cdf(ppf05,*arg) + npt.assert_almost_equal(distfn.ppf(cdf05-1e-6,*arg),ppf05, + err_msg=msg + 'ppf-cdf-median') + npt.assert_((distfn.ppf(cdf05+1e-4,*arg)>ppf05), msg + 'ppf-cdf-next') + +def check_cdf_ppf2(distfn,arg,supp,msg): + npt.assert_array_equal(distfn.ppf(distfn.cdf(supp,*arg),*arg), + supp, msg + '-roundtrip') + npt.assert_array_equal(distfn.ppf(distfn.cdf(supp,*arg)-1e-8,*arg), + supp, msg + '-roundtrip') + # -1e-8 could cause an error if pmf < 1e-8 + + +def check_cdf_ppf_private(distfn,arg,msg): + ppf05 = distfn._ppf(0.5,*arg) + cdf05 = distfn.cdf(ppf05,*arg) + npt.assert_almost_equal(distfn._ppf(cdf05-1e-6,*arg),ppf05, + err_msg=msg + '_ppf-cdf-median ') + npt.assert_((distfn._ppf(cdf05+1e-4,*arg)>ppf05), msg + '_ppf-cdf-next') + +def check_ppf_ppf(distfn, arg): + npt.assert_(distfn.ppf(0.5,*arg) < np.inf) + ppfs = distfn.ppf([0.5,0.9],*arg) + ppf_s = [distfn._ppf(0.5,*arg), distfn._ppf(0.9,*arg)] + npt.assert_(np.all(ppfs < np.inf)) + npt.assert_(ppf_s[0] == distfn.ppf(0.5,*arg)) + npt.assert_(ppf_s[1] == distfn.ppf(0.9,*arg)) + npt.assert_(ppf_s[0] == ppfs[0]) + npt.assert_(ppf_s[1] == ppfs[1]) + +def check_pmf_cdf(distfn, arg, msg): + startind = np.int(distfn._ppf(0.01,*arg)-1) + index = range(startind,startind+10) + cdfs = distfn.cdf(index,*arg) + npt.assert_almost_equal(cdfs, distfn.pmf(index, *arg).cumsum() + \ + cdfs[0] - distfn.pmf(index[0],*arg), + decimal=4, err_msg= msg + 'pmf-cdf') + +def check_generic_moment(distfn, arg, m, k, decim): + npt.assert_almost_equal(distfn.generic_moment(k,*arg), m, decimal=decim, + err_msg= str(distfn) + ' generic moment test') + +def check_moment_frozen(distfn, arg, m, k, decim): + npt.assert_almost_equal(distfn(*arg).moment(k), m, decimal=decim, + err_msg= str(distfn) + ' frozen moment test') + +def check_oth(distfn, arg, msg): + #checking other methods of distfn + meanint = round(float(distfn.stats(*arg)[0])) # closest integer to mean + npt.assert_almost_equal(distfn.sf(meanint, *arg), 1 - \ + distfn.cdf(meanint, *arg), decimal=8) + median_sf = distfn.isf(0.5, *arg) + + npt.assert_(distfn.sf(median_sf - 1, *arg) > 0.5) + npt.assert_(distfn.cdf(median_sf + 1, *arg) > 0.5) + npt.assert_equal(distfn.isf(0.5, *arg), distfn.ppf(0.5, *arg)) + +#next 3 functions copied from test_continous_extra +# adjusted + +def check_ppf_limits(distfn,arg,msg): + below,low,upp,above = distfn.ppf([-1,0,1,2], *arg) + #print distfn.name, distfn.a, low, distfn.b, upp + #print distfn.name,below,low,upp,above + assert_equal_inf_nan(distfn.a-1,low, msg + 'ppf lower bound') + assert_equal_inf_nan(distfn.b,upp, msg + 'ppf upper bound') + npt.assert_(np.isnan(below), msg + 'ppf out of bounds - below') + npt.assert_(np.isnan(above), msg + 'ppf out of bounds - above') + +def check_isf_limits(distfn,arg,msg): + below,low,upp,above = distfn.isf([-1,0,1,2], *arg) + #print distfn.name, distfn.a, low, distfn.b, upp + #print distfn.name,below,low,upp,above + assert_equal_inf_nan(distfn.a-1,upp, msg + 'isf lower bound') + assert_equal_inf_nan(distfn.b,low, msg + 'isf upper bound') + npt.assert_(np.isnan(below), msg + 'isf out of bounds - below') + npt.assert_(np.isnan(above), msg + 'isf out of bounds - above') + +def assert_equal_inf_nan(v1,v2,msg): + npt.assert_(not np.isnan(v1)) + if not np.isinf(v1): + npt.assert_almost_equal(v1, v2, decimal=10, err_msg = msg + \ + ' - finite') + else: + npt.assert_(np.isinf(v2) or np.isnan(v2), + msg + ' - infinite, v2=%s' % str(v2)) + +def check_sample_skew_kurt(distfn, arg, sk, ss, msg): + k,s = distfn.stats(moment='ks',*arg) + check_sample_meanvar, sk, k, msg + 'sample skew test' + check_sample_meanvar, ss, s, msg + 'sample kurtosis test' + + +def check_entropy(distfn,arg,msg): + ent = distfn.entropy(*arg) + #print 'Entropy =', ent + npt.assert_(not np.isnan(ent), msg + 'test Entropy is nan') + + +def check_discrete_chisquare(distfn, arg, rvs, alpha, msg): + '''perform chisquare test for random sample of a discrete distribution + + Parameters + ---------- + distname : string + name of distribution function + arg : sequence + parameters of distribution + alpha : float + significance level, threshold for p-value + + Returns + ------- + result : bool + 0 if test passes, 1 if test fails + + uses global variable debug for printing results + ''' + + # define parameters for test +## n=2000 + n = len(rvs) + nsupp = 20 + wsupp = 1.0/nsupp + +## distfn = getattr(stats, distname) +## np.random.seed(9765456) +## rvs = distfn.rvs(size=n,*arg) + + # construct intervals with minimum mass 1/nsupp + # intervalls are left-half-open as in a cdf difference + distsupport = xrange(max(distfn.a, -1000), min(distfn.b, 1000) + 1) + last = 0 + distsupp = [max(distfn.a, -1000)] + distmass = [] + for ii in distsupport: + current = distfn.cdf(ii,*arg) + if current - last >= wsupp-1e-14: + distsupp.append(ii) + distmass.append(current - last) + last = current + if current > (1-wsupp): + break + if distsupp[-1] < distfn.b: + distsupp.append(distfn.b) + distmass.append(1-last) + distsupp = np.array(distsupp) + distmass = np.array(distmass) + + # convert intervals to right-half-open as required by histogram + histsupp = distsupp+1e-8 + histsupp[0] = distfn.a + + # find sample frequencies and perform chisquare test + freq,hsupp = np.histogram(rvs,histsupp) + cdfs = distfn.cdf(distsupp,*arg) + (chis,pval) = stats.chisquare(np.array(freq),n*distmass) + + npt.assert_(pval > alpha, 'chisquare - test for %s' + ' at arg = %s with pval = %s' % (msg,str(arg),str(pval))) + + +if __name__ == "__main__": + #nose.run(argv=['', __file__]) + nose.runmodule(argv=[__file__,'-s'], exit=False) diff --git a/pywafo/src/wafo/stats/tests/test_fit.py b/pywafo/src/wafo/stats/tests/test_fit.py new file mode 100644 index 0000000..68d5793 --- /dev/null +++ b/pywafo/src/wafo/stats/tests/test_fit.py @@ -0,0 +1,98 @@ +# NOTE: contains only one test, _est_cont_fit, that is renamed so that +# nose doesn't run it +# I put this here for the record and for the case when someone wants to +# verify the quality of fit +# with current parameters: relatively small sample size, default starting values +# Ran 84 tests in 401.797s +# FAILED (failures=15) +# +#Ran 83 tests in 238.859s +#FAILED (failures=12) + +from __future__ import division, print_function, absolute_import + +import numpy.testing as npt +import numpy as np + +from wafo import stats + +from test_continuous_basic import distcont + +# this is not a proper statistical test for convergence, but only +# verifies that the estimate and true values don't differ by too much +n_repl1 = 1000 # sample size for first run +n_repl2 = 5000 # sample size for second run, if first run fails +thresh_percent = 0.25 # percent of true parameters for fail cut-off +thresh_min = 0.75 # minimum difference estimate - true to fail test + + +distslow = [ 'ncx2', 'rdist', 'gausshyper', 'recipinvgauss', 'ksone', 'genexpon', + 'vonmises', 'rice', 'mielke', + 'powerlognorm', 'kstwobign', 'tukeylambda','betaprime', 'gengamma', + 'johnsonsb', 'burr', 'truncexpon', 'pearson3', 'exponweib', 'nakagami', + 'wrapcauchy'] +dist_rarely_fitted = ['f', 'ncf', 'nct', 'chi'] +distskip = distslow + dist_rarely_fitted + +#distcont = [['genextreme', (3.3184017469423535,)]] +#@npt.dec.slow +def test_cont_fit(): + # this tests the closeness of the estimated parameters to the true + # parameters with fit method of continuous distributions + for distname, arg in distcont: + if distname not in distskip: + yield check_cont_fit, distname,arg + +@npt.dec.slow +def _est_cont_fit_slow(): + # this tests the closeness of the estimated parameters to the true + # parameters with fit method of continuous distributions + # Note: is slow, some distributions don't converge with sample size <= 10000 + for distname, arg in distcont: + if distname in distslow: + yield check_cont_fit, distname,arg + +def test_lognorm_fit_ticket1131(): + params = [(2.1, 1.,1.), (1.,10.,1.), (1.,1.,10.)] + for param in params: + yield check_cont_fit, 'lognorm', param + +def check_cont_fit(distname,arg): + distfn = getattr(stats, distname) + rvs = distfn.rvs(size=n_repl1,*arg) + est = distfn.fit(rvs) #, *arg) # start with default values + n = distfn.numargs + 2 + truearg = np.hstack([arg,[0.0, 1.0]])[:n] + + diff = est-truearg + + txt = '' + diffthreshold = np.max(np.vstack([truearg*thresh_percent, + np.ones(distfn.numargs+2)*thresh_min]),0) + # threshold for location + diffthreshold[-2] = np.max([np.abs(rvs.mean())*thresh_percent,thresh_min]) + + if np.any(np.isnan(est)): + raise AssertionError('nan returned in fit') + else: + if np.any((np.abs(diff) - diffthreshold) > 0.0): +## txt = 'WARNING - diff too large with small sample' +## print 'parameter diff =', diff - diffthreshold, txt + rvs = np.concatenate([rvs,distfn.rvs(size=n_repl2-n_repl1,*arg)]) + est = distfn.fit(rvs) #,*arg) + truearg = np.hstack([arg,[0.0,1.0]])[:n] + diff = est-truearg + if np.any((np.abs(diff) - diffthreshold) > 0.0): + txt = 'parameter: %s\n' % str(truearg) + txt += 'estimated: %s\n' % str(est) + txt += 'diff : %s\n' % str(diff) + raise AssertionError('fit not very good in %s\n' % distfn.name + txt) + + +if __name__ == "__main__": + check_cont_fit('bradford', (0.29891359763170633,)) +# check_cont_fit('lognorm', (10,1,1)) +# check_cont_fit('ncx2', (21, 1.0560465975116415)) + import nose + #nose.run(argv=['', __file__]) + nose.runmodule(argv=[__file__,'-s'], exit=False) diff --git a/pywafo/src/wafo/stats/tests/test_kdeoth.py b/pywafo/src/wafo/stats/tests/test_kdeoth.py new file mode 100644 index 0000000..ee822b7 --- /dev/null +++ b/pywafo/src/wafo/stats/tests/test_kdeoth.py @@ -0,0 +1,36 @@ + + + +from wafo import stats +import numpy as np +from numpy.testing import assert_almost_equal, assert_ + +def test_kde_1d(): + #some basic tests comparing to normal distribution + np.random.seed(8765678) + n_basesample = 500 + xn = np.random.randn(n_basesample) + xnmean = xn.mean() + xnstd = xn.std(ddof=1) + + # get kde for original sample + gkde = stats.gaussian_kde(xn) + + # evaluate the density funtion for the kde for some points + xs = np.linspace(-7,7,501) + kdepdf = gkde.evaluate(xs) + normpdf = stats.norm.pdf(xs, loc=xnmean, scale=xnstd) + intervall = xs[1] - xs[0] + + assert_(np.sum((kdepdf - normpdf)**2)*intervall < 0.01) + prob1 = gkde.integrate_box_1d(xnmean, np.inf) + prob2 = gkde.integrate_box_1d(-np.inf, xnmean) + assert_almost_equal(prob1, 0.5, decimal=1) + assert_almost_equal(prob2, 0.5, decimal=1) + assert_almost_equal(gkde.integrate_box(xnmean, np.inf), prob1, decimal=13) + assert_almost_equal(gkde.integrate_box(-np.inf, xnmean), prob2, decimal=13) + + assert_almost_equal(gkde.integrate_kde(gkde), + (kdepdf**2).sum()*intervall, decimal=2) + assert_almost_equal(gkde.integrate_gaussian(xnmean, xnstd**2), + (kdepdf*normpdf).sum()*intervall, decimal=2) diff --git a/pywafo/src/wafo/stats/tests/test_mstats_basic.py b/pywafo/src/wafo/stats/tests/test_mstats_basic.py new file mode 100644 index 0000000..ca45e8b --- /dev/null +++ b/pywafo/src/wafo/stats/tests/test_mstats_basic.py @@ -0,0 +1,490 @@ +""" +Tests for the stats.mstats module (support for maskd arrays) +""" + + +import numpy as np +from numpy import nan +import numpy.ma as ma +from numpy.ma import masked, nomask + +import scipy.stats.mstats as mstats +from numpy.testing import TestCase, run_module_suite +from numpy.ma.testutils import assert_equal, assert_almost_equal, \ + assert_array_almost_equal, assert_ + + +class TestMquantiles(TestCase): + """Regression tests for mstats module.""" + def test_mquantiles_limit_keyword(self): + """Ticket #867""" + data = np.array([[ 6., 7., 1.], + [ 47., 15., 2.], + [ 49., 36., 3.], + [ 15., 39., 4.], + [ 42., 40., -999.], + [ 41., 41., -999.], + [ 7., -999., -999.], + [ 39., -999., -999.], + [ 43., -999., -999.], + [ 40., -999., -999.], + [ 36., -999., -999.]]) + desired = [[19.2, 14.6, 1.45], + [40.0, 37.5, 2.5 ], + [42.8, 40.05, 3.55]] + quants = mstats.mquantiles(data, axis=0, limit=(0, 50)) + assert_almost_equal(quants, desired) + + + +class TestGMean(TestCase): + def test_1D(self): + a = (1,2,3,4) + actual= mstats.gmean(a) + desired = np.power(1*2*3*4,1./4.) + assert_almost_equal(actual, desired,decimal=14) + + desired1 = mstats.gmean(a,axis=-1) + assert_almost_equal(actual, desired1, decimal=14) + assert_(not isinstance(desired1, ma.MaskedArray)) + # + a = ma.array((1,2,3,4),mask=(0,0,0,1)) + actual= mstats.gmean(a) + desired = np.power(1*2*3,1./3.) + assert_almost_equal(actual, desired,decimal=14) + + desired1 = mstats.gmean(a,axis=-1) + assert_almost_equal(actual, desired1, decimal=14) + # + def test_2D(self): + a = ma.array(((1,2,3,4),(1,2,3,4),(1,2,3,4)), + mask=((0,0,0,0),(1,0,0,1),(0,1,1,0))) + actual= mstats.gmean(a) + desired = np.array((1,2,3,4)) + assert_array_almost_equal(actual, desired, decimal=14) + # + desired1 = mstats.gmean(a,axis=0) + assert_array_almost_equal(actual, desired1, decimal=14) + # + actual= mstats.gmean(a, -1) + desired = ma.array((np.power(1*2*3*4,1./4.), + np.power(2*3,1./2.), + np.power(1*4,1./2.))) + assert_array_almost_equal(actual, desired, decimal=14) + +class TestHMean(TestCase): + def test_1D(self): + a = (1,2,3,4) + actual= mstats.hmean(a) + desired = 4. / (1./1 + 1./2 + 1./3 + 1./4) + assert_almost_equal(actual, desired, decimal=14) + desired1 = mstats.hmean(ma.array(a),axis=-1) + assert_almost_equal(actual, desired1, decimal=14) + # + a = ma.array((1,2,3,4),mask=(0,0,0,1)) + actual= mstats.hmean(a) + desired = 3. / (1./1 + 1./2 + 1./3) + assert_almost_equal(actual, desired,decimal=14) + desired1 = mstats.hmean(a,axis=-1) + assert_almost_equal(actual, desired1, decimal=14) + + def test_2D(self): + a = ma.array(((1,2,3,4),(1,2,3,4),(1,2,3,4)), + mask=((0,0,0,0),(1,0,0,1),(0,1,1,0))) + actual= mstats.hmean(a) + desired = ma.array((1,2,3,4)) + assert_array_almost_equal(actual, desired, decimal=14) + # + actual1 = mstats.hmean(a,axis=-1) + desired = (4./(1/1.+1/2.+1/3.+1/4.), + 2./(1/2.+1/3.), + 2./(1/1.+1/4.) + ) + assert_array_almost_equal(actual1, desired, decimal=14) + + +class TestRanking(TestCase): + # + def __init__(self, *args, **kwargs): + TestCase.__init__(self, *args, **kwargs) + # + def test_ranking(self): + x = ma.array([0,1,1,1,2,3,4,5,5,6,]) + assert_almost_equal(mstats.rankdata(x),[1,3,3,3,5,6,7,8.5,8.5,10]) + x[[3,4]] = masked + assert_almost_equal(mstats.rankdata(x),[1,2.5,2.5,0,0,4,5,6.5,6.5,8]) + assert_almost_equal(mstats.rankdata(x,use_missing=True), + [1,2.5,2.5,4.5,4.5,4,5,6.5,6.5,8]) + x = ma.array([0,1,5,1,2,4,3,5,1,6,]) + assert_almost_equal(mstats.rankdata(x),[1,3,8.5,3,5,7,6,8.5,3,10]) + x = ma.array([[0,1,1,1,2], [3,4,5,5,6,]]) + assert_almost_equal(mstats.rankdata(x),[[1,3,3,3,5],[6,7,8.5,8.5,10]]) + assert_almost_equal(mstats.rankdata(x,axis=1),[[1,3,3,3,5],[1,2,3.5,3.5,5]]) + assert_almost_equal(mstats.rankdata(x,axis=0),[[1,1,1,1,1],[2,2,2,2,2,]]) + + +class TestCorr(TestCase): + # + def test_pearsonr(self): + "Tests some computations of Pearson's r" + x = ma.arange(10) + olderr = np.seterr(all='ignore') + try: + assert_almost_equal(mstats.pearsonr(x,x)[0], 1.0) + assert_almost_equal(mstats.pearsonr(x,x[::-1])[0], -1.0) + + x = ma.array(x, mask=True) + pr = mstats.pearsonr(x,x) + finally: + np.seterr(**olderr) + assert_(pr[0] is masked) + assert_(pr[1] is masked) + # + def test_spearmanr(self): + "Tests some computations of Spearman's rho" + (x, y) = ([5.05,6.75,3.21,2.66],[1.65,2.64,2.64,6.95]) + assert_almost_equal(mstats.spearmanr(x,y)[0], -0.6324555) + (x, y) = ([5.05,6.75,3.21,2.66,np.nan],[1.65,2.64,2.64,6.95,np.nan]) + (x, y) = (ma.fix_invalid(x), ma.fix_invalid(y)) + assert_almost_equal(mstats.spearmanr(x,y)[0], -0.6324555) + # + x = [ 2.0, 47.4, 42.0, 10.8, 60.1, 1.7, 64.0, 63.1, + 1.0, 1.4, 7.9, 0.3, 3.9, 0.3, 6.7] + y = [22.6, 08.3, 44.4, 11.9, 24.6, 0.6, 5.7, 41.6, + 0.0, 0.6, 6.7, 3.8, 1.0, 1.2, 1.4] + assert_almost_equal(mstats.spearmanr(x,y)[0], 0.6887299) + x = [ 2.0, 47.4, 42.0, 10.8, 60.1, 1.7, 64.0, 63.1, + 1.0, 1.4, 7.9, 0.3, 3.9, 0.3, 6.7, np.nan] + y = [22.6, 08.3, 44.4, 11.9, 24.6, 0.6, 5.7, 41.6, + 0.0, 0.6, 6.7, 3.8, 1.0, 1.2, 1.4, np.nan] + (x, y) = (ma.fix_invalid(x), ma.fix_invalid(y)) + assert_almost_equal(mstats.spearmanr(x,y)[0], 0.6887299) + # + def test_kendalltau(self): + "Tests some computations of Kendall's tau" + x = ma.fix_invalid([5.05, 6.75, 3.21, 2.66,np.nan]) + y = ma.fix_invalid([1.65, 26.5, -5.93, 7.96, np.nan]) + z = ma.fix_invalid([1.65, 2.64, 2.64, 6.95, np.nan]) + assert_almost_equal(np.asarray(mstats.kendalltau(x,y)), + [+0.3333333,0.4969059]) + assert_almost_equal(np.asarray(mstats.kendalltau(x,z)), + [-0.5477226,0.2785987]) + # + x = ma.fix_invalid([ 0, 0, 0, 0,20,20, 0,60, 0,20, + 10,10, 0,40, 0,20, 0, 0, 0, 0, 0, np.nan]) + y = ma.fix_invalid([ 0,80,80,80,10,33,60, 0,67,27, + 25,80,80,80,80,80,80, 0,10,45, np.nan, 0]) + result = mstats.kendalltau(x,y) + assert_almost_equal(np.asarray(result), [-0.1585188, 0.4128009]) + # + def test_kendalltau_seasonal(self): + "Tests the seasonal Kendall tau." + x = [[nan,nan, 4, 2, 16, 26, 5, 1, 5, 1, 2, 3, 1], + [ 4, 3, 5, 3, 2, 7, 3, 1, 1, 2, 3, 5, 3], + [ 3, 2, 5, 6, 18, 4, 9, 1, 1,nan, 1, 1,nan], + [nan, 6, 11, 4, 17,nan, 6, 1, 1, 2, 5, 1, 1]] + x = ma.fix_invalid(x).T + output = mstats.kendalltau_seasonal(x) + assert_almost_equal(output['global p-value (indep)'], 0.008, 3) + assert_almost_equal(output['seasonal p-value'].round(2), + [0.18,0.53,0.20,0.04]) + # + def test_pointbiserial(self): + "Tests point biserial" + x = [1,0,1,1,1,1,0,1,0,0,0,1,1,0,0,0,1,1,1,0,0,0,0,0,0,0,0,1,0, + 0,0,0,0,1,-1] + y = [14.8,13.8,12.4,10.1,7.1,6.1,5.8,4.6,4.3,3.5,3.3,3.2,3.0, + 2.8,2.8,2.5,2.4,2.3,2.1,1.7,1.7,1.5,1.3,1.3,1.2,1.2,1.1, + 0.8,0.7,0.6,0.5,0.2,0.2,0.1,np.nan] + assert_almost_equal(mstats.pointbiserialr(x, y)[0], 0.36149, 5) + + +class TestTrimming(TestCase): + # + def test_trim(self): + "Tests trimming" + a = ma.arange(10) + assert_equal(mstats.trim(a), [0,1,2,3,4,5,6,7,8,9]) + a = ma.arange(10) + assert_equal(mstats.trim(a,(2,8)), [None,None,2,3,4,5,6,7,8,None]) + a = ma.arange(10) + assert_equal(mstats.trim(a,limits=(2,8),inclusive=(False,False)), + [None,None,None,3,4,5,6,7,None,None]) + a = ma.arange(10) + assert_equal(mstats.trim(a,limits=(0.1,0.2),relative=True), + [None,1,2,3,4,5,6,7,None,None]) + # + a = ma.arange(12) + a[[0,-1]] = a[5] = masked + assert_equal(mstats.trim(a,(2,8)), + [None,None,2,3,4,None,6,7,8,None,None,None]) + # + x = ma.arange(100).reshape(10,10) + trimx = mstats.trim(x,(0.1,0.2),relative=True,axis=None) + assert_equal(trimx._mask.ravel(),[1]*10+[0]*70+[1]*20) + trimx = mstats.trim(x,(0.1,0.2),relative=True,axis=0) + assert_equal(trimx._mask.ravel(),[1]*10+[0]*70+[1]*20) + trimx = mstats.trim(x,(0.1,0.2),relative=True,axis=-1) + assert_equal(trimx._mask.T.ravel(),[1]*10+[0]*70+[1]*20) + # + x = ma.arange(110).reshape(11,10) + x[1] = masked + trimx = mstats.trim(x,(0.1,0.2),relative=True,axis=None) + assert_equal(trimx._mask.ravel(),[1]*20+[0]*70+[1]*20) + trimx = mstats.trim(x,(0.1,0.2),relative=True,axis=0) + assert_equal(trimx._mask.ravel(),[1]*20+[0]*70+[1]*20) + trimx = mstats.trim(x.T,(0.1,0.2),relative=True,axis=-1) + assert_equal(trimx.T._mask.ravel(),[1]*20+[0]*70+[1]*20) + # + def test_trim_old(self): + "Tests trimming." + x = ma.arange(100) + assert_equal(mstats.trimboth(x).count(), 60) + assert_equal(mstats.trimtail(x,tail='r').count(), 80) + x[50:70] = masked + trimx = mstats.trimboth(x) + assert_equal(trimx.count(), 48) + assert_equal(trimx._mask, [1]*16 + [0]*34 + [1]*20 + [0]*14 + [1]*16) + x._mask = nomask + x.shape = (10,10) + assert_equal(mstats.trimboth(x).count(), 60) + assert_equal(mstats.trimtail(x).count(), 80) + # + def test_trimmedmean(self): + "Tests the trimmed mean." + data = ma.array([ 77, 87, 88,114,151,210,219,246,253,262, + 296,299,306,376,428,515,666,1310,2611]) + assert_almost_equal(mstats.trimmed_mean(data,0.1), 343, 0) + assert_almost_equal(mstats.trimmed_mean(data,(0.1,0.1)), 343, 0) + assert_almost_equal(mstats.trimmed_mean(data,(0.2,0.2)), 283, 0) + # + def test_trimmed_stde(self): + "Tests the trimmed mean standard error." + data = ma.array([ 77, 87, 88,114,151,210,219,246,253,262, + 296,299,306,376,428,515,666,1310,2611]) + assert_almost_equal(mstats.trimmed_stde(data,(0.2,0.2)), 56.13193, 5) + assert_almost_equal(mstats.trimmed_stde(data,0.2), 56.13193, 5) + # + def test_winsorization(self): + "Tests the Winsorization of the data." + data = ma.array([ 77, 87, 88,114,151,210,219,246,253,262, + 296,299,306,376,428,515,666,1310,2611]) + assert_almost_equal(mstats.winsorize(data,(0.2,0.2)).var(ddof=1), + 21551.4, 1) + data[5] = masked + winsorized = mstats.winsorize(data) + assert_equal(winsorized.mask, data.mask) + + +class TestMoments(TestCase): + """ + Comparison numbers are found using R v.1.5.1 + note that length(testcase) = 4 + testmathworks comes from documentation for the + Statistics Toolbox for Matlab and can be found at both + http://www.mathworks.com/access/helpdesk/help/toolbox/stats/kurtosis.shtml + http://www.mathworks.com/access/helpdesk/help/toolbox/stats/skewness.shtml + Note that both test cases came from here. + """ + testcase = [1,2,3,4] + testmathworks = ma.fix_invalid([1.165 , 0.6268, 0.0751, 0.3516, -0.6965, + np.nan]) + def test_moment(self): + """ + mean((testcase-mean(testcase))**power,axis=0),axis=0))**power))""" + y = mstats.moment(self.testcase,1) + assert_almost_equal(y,0.0,10) + y = mstats.moment(self.testcase,2) + assert_almost_equal(y,1.25) + y = mstats.moment(self.testcase,3) + assert_almost_equal(y,0.0) + y = mstats.moment(self.testcase,4) + assert_almost_equal(y,2.5625) + def test_variation(self): + """variation = samplestd/mean """ +## y = stats.variation(self.shoes[0]) +## assert_almost_equal(y,21.8770668) + y = mstats.variation(self.testcase) + assert_almost_equal(y,0.44721359549996, 10) + + def test_skewness(self): + """ + sum((testmathworks-mean(testmathworks,axis=0))**3,axis=0)/((sqrt(var(testmathworks)*4/5))**3)/5 + """ + y = mstats.skew(self.testmathworks) + assert_almost_equal(y,-0.29322304336607,10) + y = mstats.skew(self.testmathworks,bias=0) + assert_almost_equal(y,-0.437111105023940,10) + y = mstats.skew(self.testcase) + assert_almost_equal(y,0.0,10) + + def test_kurtosis(self): + """ + sum((testcase-mean(testcase,axis=0))**4,axis=0)/((sqrt(var(testcase)*3/4))**4)/4 + sum((test2-mean(testmathworks,axis=0))**4,axis=0)/((sqrt(var(testmathworks)*4/5))**4)/5 + Set flags for axis = 0 and + fisher=0 (Pearson's definition of kurtosis for compatibility with Matlab) + """ + y = mstats.kurtosis(self.testmathworks,0,fisher=0,bias=1) + assert_almost_equal(y, 2.1658856802973,10) + # Note that MATLAB has confusing docs for the following case + # kurtosis(x,0) gives an unbiased estimate of Pearson's skewness + # kurtosis(x) gives a biased estimate of Fisher's skewness (Pearson-3) + # The MATLAB docs imply that both should give Fisher's + y = mstats.kurtosis(self.testmathworks,fisher=0,bias=0) + assert_almost_equal(y, 3.663542721189047,10) + y = mstats.kurtosis(self.testcase,0,0) + assert_almost_equal(y,1.64) + # + def test_mode(self): + "Tests the mode" + # + a1 = [0,0,0,1,1,1,2,3,3,3,3,4,5,6,7] + a2 = np.reshape(a1, (3,5)) + ma1 = ma.masked_where(ma.array(a1)>2,a1) + ma2 = ma.masked_where(a2>2, a2) + assert_equal(mstats.mode(a1, axis=None), (3,4)) + assert_equal(mstats.mode(ma1, axis=None), (0,3)) + assert_equal(mstats.mode(a2, axis=None), (3,4)) + assert_equal(mstats.mode(ma2, axis=None), (0,3)) + assert_equal(mstats.mode(a2, axis=0), ([[0,0,0,1,1]],[[1,1,1,1,1]])) + assert_equal(mstats.mode(ma2, axis=0), ([[0,0,0,1,1]],[[1,1,1,1,1]])) + assert_equal(mstats.mode(a2, axis=-1), ([[0],[3],[3]], [[3],[3],[1]])) + assert_equal(mstats.mode(ma2, axis=-1), ([[0],[1],[0]], [[3],[1],[0]])) + + +class TestPercentile(TestCase): + def setUp(self): + self.a1 = [3,4,5,10,-3,-5,6] + self.a2 = [3,-6,-2,8,7,4,2,1] + self.a3 = [3.,4,5,10,-3,-5,-6,7.0] + + def test_percentile(self): + x = np.arange(8) * 0.5 + assert_equal(mstats.scoreatpercentile(x, 0), 0.) + assert_equal(mstats.scoreatpercentile(x, 100), 3.5) + assert_equal(mstats.scoreatpercentile(x, 50), 1.75) + + def test_2D(self): + x = ma.array([[1, 1, 1], + [1, 1, 1], + [4, 4, 3], + [1, 1, 1], + [1, 1, 1]]) + assert_equal(mstats.scoreatpercentile(x,50), [1,1,1]) + + +class TestVariability(TestCase): + """ Comparison numbers are found using R v.1.5.1 + note that length(testcase) = 4 + """ + testcase = ma.fix_invalid([1,2,3,4,np.nan]) + + def test_signaltonoise(self): + """ + this is not in R, so used + mean(testcase,axis=0)/(sqrt(var(testcase)*3/4)) """ + #y = stats.signaltonoise(self.shoes[0]) + #assert_approx_equal(y,4.5709967) + y = mstats.signaltonoise(self.testcase) + assert_almost_equal(y,2.236067977) + + def test_sem(self): + """ + this is not in R, so used + sqrt(var(testcase)*3/4)/sqrt(3) + """ + #y = stats.sem(self.shoes[0]) + #assert_approx_equal(y,0.775177399) + y = mstats.sem(self.testcase) + assert_almost_equal(y,0.6454972244) + + def test_zmap(self): + """ + not in R, so tested by using + (testcase[i]-mean(testcase,axis=0))/sqrt(var(testcase)*3/4) + """ + y = mstats.zmap(self.testcase, self.testcase) + desired_unmaskedvals = ([-1.3416407864999, -0.44721359549996 , + 0.44721359549996 , 1.3416407864999]) + assert_array_almost_equal(desired_unmaskedvals, + y.data[y.mask==False], decimal=12) + + def test_zscore(self): + """ + not in R, so tested by using + (testcase[i]-mean(testcase,axis=0))/sqrt(var(testcase)*3/4) + """ + y = mstats.zscore(self.testcase) + desired = ma.fix_invalid([-1.3416407864999, -0.44721359549996 , + 0.44721359549996 , 1.3416407864999, np.nan]) + assert_almost_equal(desired, y, decimal=12) + + + +class TestMisc(TestCase): + # + def test_obrientransform(self): + "Tests Obrien transform" + args = [[5]*5+[6]*11+[7]*9+[8]*3+[9]*2+[10]*2, + [6]+[7]*2+[8]*4+[9]*9+[10]*16] + result = [5*[3.1828]+11*[0.5591]+9*[0.0344]+3*[1.6086]+2*[5.2817]+2*[11.0538], + [10.4352]+2*[4.8599]+4*[1.3836]+9*[0.0061]+16*[0.7277]] + assert_almost_equal(np.round(mstats.obrientransform(*args).T,4), + result,4) + # + def test_kstwosamp(self): + "Tests the Kolmogorov-Smirnov 2 samples test" + x = [[nan,nan, 4, 2, 16, 26, 5, 1, 5, 1, 2, 3, 1], + [ 4, 3, 5, 3, 2, 7, 3, 1, 1, 2, 3, 5, 3], + [ 3, 2, 5, 6, 18, 4, 9, 1, 1,nan, 1, 1,nan], + [nan, 6, 11, 4, 17,nan, 6, 1, 1, 2, 5, 1, 1]] + x = ma.fix_invalid(x).T + (winter,spring,summer,fall) = x.T + # + assert_almost_equal(np.round(mstats.ks_twosamp(winter,spring),4), + (0.1818,0.9892)) + assert_almost_equal(np.round(mstats.ks_twosamp(winter,spring,'g'),4), + (0.1469,0.7734)) + assert_almost_equal(np.round(mstats.ks_twosamp(winter,spring,'l'),4), + (0.1818,0.6744)) + # + def test_friedmanchisq(self): + "Tests the Friedman Chi-square test" + # No missing values + args = ([9.0,9.5,5.0,7.5,9.5,7.5,8.0,7.0,8.5,6.0], + [7.0,6.5,7.0,7.5,5.0,8.0,6.0,6.5,7.0,7.0], + [6.0,8.0,4.0,6.0,7.0,6.5,6.0,4.0,6.5,3.0]) + result = mstats.friedmanchisquare(*args) + assert_almost_equal(result[0], 10.4737, 4) + assert_almost_equal(result[1], 0.005317, 6) + # Missing values + x = [[nan,nan, 4, 2, 16, 26, 5, 1, 5, 1, 2, 3, 1], + [ 4, 3, 5, 3, 2, 7, 3, 1, 1, 2, 3, 5, 3], + [ 3, 2, 5, 6, 18, 4, 9, 1, 1,nan, 1, 1,nan], + [nan, 6, 11, 4, 17,nan, 6, 1, 1, 2, 5, 1, 1]] + x = ma.fix_invalid(x) + result = mstats.friedmanchisquare(*x) + assert_almost_equal(result[0], 2.0156, 4) + assert_almost_equal(result[1], 0.5692, 4) + + +def test_regress_simple(): + """Regress a line with sinusoidal noise. Test for #1273.""" + x = np.linspace(0, 100, 100) + y = 0.2 * np.linspace(0, 100, 100) + 10 + y += np.sin(np.linspace(0, 20, 100)) + + slope, intercept, r_value, p_value, sterr = mstats.linregress(x, y) + assert_almost_equal(slope, 0.19644990055858422) + assert_almost_equal(intercept, 10.211269918932341) + + +def test_plotting_positions(): + """Regression test for #1256""" + pos = mstats.plotting_positions(np.arange(3), 0, 0) + assert_array_almost_equal(pos.data, np.array([0.25, 0.5, 0.75])) + + +if __name__ == "__main__": + run_module_suite() diff --git a/pywafo/src/wafo/stats/tests/test_mstats_extras.py b/pywafo/src/wafo/stats/tests/test_mstats_extras.py new file mode 100644 index 0000000..2fe5d72 --- /dev/null +++ b/pywafo/src/wafo/stats/tests/test_mstats_extras.py @@ -0,0 +1,103 @@ +# pylint: disable-msg=W0611, W0612, W0511,R0201 +"""Tests suite for maskedArray statistics. + +:author: Pierre Gerard-Marchant +:contact: pierregm_at_uga_dot_edu +""" +__author__ = "Pierre GF Gerard-Marchant ($Author: backtopop $)" + +import numpy as np + +import numpy.ma as ma + +import scipy.stats.mstats as ms +#import scipy.stats.mmorestats as mms + +from numpy.testing import TestCase, run_module_suite, assert_equal, \ + assert_almost_equal, assert_ + + +class TestMisc(TestCase): + # + def __init__(self, *args, **kwargs): + TestCase.__init__(self, *args, **kwargs) + # + def test_mjci(self): + "Tests the Marits-Jarrett estimator" + data = ma.array([ 77, 87, 88,114,151,210,219,246,253,262, + 296,299,306,376,428,515,666,1310,2611]) + assert_almost_equal(ms.mjci(data),[55.76819,45.84028,198.87875],5) + # + def test_trimmedmeanci(self): + "Tests the confidence intervals of the trimmed mean." + data = ma.array([545,555,558,572,575,576,578,580, + 594,605,635,651,653,661,666]) + assert_almost_equal(ms.trimmed_mean(data,0.2), 596.2, 1) + assert_equal(np.round(ms.trimmed_mean_ci(data,(0.2,0.2)),1), + [561.8, 630.6]) + # + def test_idealfourths(self): + "Tests ideal-fourths" + test = np.arange(100) + assert_almost_equal(np.asarray(ms.idealfourths(test)), + [24.416667,74.583333],6) + test_2D = test.repeat(3).reshape(-1,3) + assert_almost_equal(ms.idealfourths(test_2D, axis=0), + [[24.416667,24.416667,24.416667], + [74.583333,74.583333,74.583333]],6) + assert_almost_equal(ms.idealfourths(test_2D, axis=1), + test.repeat(2).reshape(-1,2)) + test = [0,0] + _result = ms.idealfourths(test) + assert_(np.isnan(_result).all()) + +#.............................................................................. +class TestQuantiles(TestCase): + # + def __init__(self, *args, **kwargs): + TestCase.__init__(self, *args, **kwargs) + # + def test_hdquantiles(self): + data = [0.706560797,0.727229578,0.990399276,0.927065621,0.158953014, + 0.887764025,0.239407086,0.349638551,0.972791145,0.149789972, + 0.936947700,0.132359948,0.046041972,0.641675031,0.945530547, + 0.224218684,0.771450991,0.820257774,0.336458052,0.589113496, + 0.509736129,0.696838829,0.491323573,0.622767425,0.775189248, + 0.641461450,0.118455200,0.773029450,0.319280007,0.752229111, + 0.047841438,0.466295911,0.583850781,0.840581845,0.550086491, + 0.466470062,0.504765074,0.226855960,0.362641207,0.891620942, + 0.127898691,0.490094097,0.044882048,0.041441695,0.317976349, + 0.504135618,0.567353033,0.434617473,0.636243375,0.231803616, + 0.230154113,0.160011327,0.819464108,0.854706985,0.438809221, + 0.487427267,0.786907310,0.408367937,0.405534192,0.250444460, + 0.995309248,0.144389588,0.739947527,0.953543606,0.680051621, + 0.388382017,0.863530727,0.006514031,0.118007779,0.924024803, + 0.384236354,0.893687694,0.626534881,0.473051932,0.750134705, + 0.241843555,0.432947602,0.689538104,0.136934797,0.150206859, + 0.474335206,0.907775349,0.525869295,0.189184225,0.854284286, + 0.831089744,0.251637345,0.587038213,0.254475554,0.237781276, + 0.827928620,0.480283781,0.594514455,0.213641488,0.024194386, + 0.536668589,0.699497811,0.892804071,0.093835427,0.731107772] + # + assert_almost_equal(ms.hdquantiles(data,[0., 1.]), + [0.006514031, 0.995309248]) + hdq = ms.hdquantiles(data,[0.25, 0.5, 0.75]) + assert_almost_equal(hdq, [0.253210762, 0.512847491, 0.762232442,]) + hdq = ms.hdquantiles_sd(data,[0.25, 0.5, 0.75]) + assert_almost_equal(hdq, [0.03786954, 0.03805389, 0.03800152,], 4) + # + data = np.array(data).reshape(10,10) + hdq = ms.hdquantiles(data,[0.25,0.5,0.75],axis=0) + assert_almost_equal(hdq[:,0], ms.hdquantiles(data[:,0],[0.25,0.5,0.75])) + assert_almost_equal(hdq[:,-1], ms.hdquantiles(data[:,-1],[0.25,0.5,0.75])) + hdq = ms.hdquantiles(data,[0.25,0.5,0.75],axis=0,var=True) + assert_almost_equal(hdq[...,0], + ms.hdquantiles(data[:,0],[0.25,0.5,0.75],var=True)) + assert_almost_equal(hdq[...,-1], + ms.hdquantiles(data[:,-1],[0.25,0.5,0.75], var=True)) + + +############################################################################### + +if __name__ == "__main__": + run_module_suite() diff --git a/pywafo/src/wafo/stats/tests/test_stats.py b/pywafo/src/wafo/stats/tests/test_stats.py new file mode 100644 index 0000000..e716d7d --- /dev/null +++ b/pywafo/src/wafo/stats/tests/test_stats.py @@ -0,0 +1,2205 @@ +""" Test functions for stats module + + WRITTEN BY LOUIS LUANGKESORN FOR THE STATS MODULE + BASED ON WILKINSON'S STATISTICS QUIZ + http://www.stanford.edu/~clint/bench/wilk.txt + + Additional tests by a host of SciPy developers. +""" +from __future__ import division, print_function, absolute_import + +import sys + +from numpy.testing import TestCase, assert_, assert_equal, \ + assert_almost_equal, assert_array_almost_equal, assert_array_equal, \ + assert_approx_equal, assert_raises, run_module_suite, \ + assert_allclose, dec +from numpy import array, arange, float32, float64, power +import numpy as np + +import wafo.stats as stats + + +""" Numbers in docstrings begining with 'W' refer to the section numbers + and headings found in the STATISTICS QUIZ of Leland Wilkinson. These are + considered to be essential functionality. True testing and + evaluation of a statistics package requires use of the + NIST Statistical test data. See McCoullough(1999) Assessing The Reliability + of Statistical Software for a test methodology and its + implementation in testing SAS, SPSS, and S-Plus +""" + +## Datasets +## These data sets are from the nasty.dat sets used by Wilkinson +## for MISS, need to be able to represent missing values +## For completeness, I should write the relevant tests and count them as failures +## Somewhat acceptable, since this is still beta software. It would count as a +## good target for 1.0 status +X = array([1,2,3,4,5,6,7,8,9],float) +ZERO= array([0,0,0,0,0,0,0,0,0], float) +#MISS=array([.,.,.,.,.,.,.,.,.], float) +BIG=array([99999991,99999992,99999993,99999994,99999995,99999996,99999997,99999998,99999999],float) +LITTLE=array([0.99999991,0.99999992,0.99999993,0.99999994,0.99999995,0.99999996,0.99999997,0.99999998,0.99999999],float) +HUGE=array([1e+12,2e+12,3e+12,4e+12,5e+12,6e+12,7e+12,8e+12,9e+12],float) +TINY=array([1e-12,2e-12,3e-12,4e-12,5e-12,6e-12,7e-12,8e-12,9e-12],float) +ROUND=array([0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5],float) +X2 = X * X +X3 = X2 * X +X4 = X3 * X +X5 = X4 * X +X6 = X5 * X +X7 = X6 * X +X8 = X7 * X +X9 = X8 * X + +class TestRound(TestCase): + """ W.II. ROUND + + You should get the numbers 1 to 9. Many language compilers, + such as Turbo Pascal and Lattice C, fail this test (they round + numbers inconsistently). Needless to say, statical packages + written in these languages may fail the test as well. You can + also check the following expressions: + Y = INT(2.6*7 -0.2) (Y should be 18) + Y = 2-INT(EXP(LOG(SQR(2)*SQR(2)))) (Y should be 0) + Y = INT(3-EXP(LOG(SQR(2)*SQR(2)))) (Y should be 1) + INT is the integer function. It converts decimal numbers to + integers by throwing away numbers after the decimal point. EXP + is exponential, LOG is logarithm, and SQR is suqare root. You may + have to substitute similar names for these functions for different + packages. Since the square of a square root should return the same + number, and the exponential of a log should return the same number, + we should get back a 2 from this function of functions. By taking + the integer result and subtracting from 2, we are exposing the + roundoff errors. These simple functions are at the heart of + statistical calculations. + """ + + def test_rounding0(self): + """ W.II.A.0. Print ROUND with only one digit. + + You should get the numbers 1 to 9. Many language compilers, + such as Turbo Pascal and Lattice C, fail this test (they round + numbers inconsistently). Needless to say, statical packages + written in these languages may fail the test as well. + """ + if sys.version_info[0] >= 3: + # round to even + for i in range(0,9): + y = round(ROUND[i]) + assert_equal(y, 2*((i+1)//2)) + else: + for i in range(0,9): + y = round(ROUND[i]) + assert_equal(y,i+1) + + def test_rounding1(self): + """ W.II.A.1. Y = INT(2.6*7 -0.2) (Y should be 18)""" + y = int(2.6*7 -0.2) + assert_equal(y, 18) + + def test_rounding2(self): + """ W.II.A.2. Y = 2-INT(EXP(LOG(SQR(2)*SQR(2)))) (Y should be 0)""" + y=2-int(np.exp(np.log(np.sqrt(2.)*np.sqrt(2.)))) + assert_equal(y,0) + + def test_rounding3(self): + """ W.II.A.3. Y = INT(3-EXP(LOG(SQR(2)*SQR(2)))) (Y should be 1)""" + y=(int(round((3-np.exp(np.log(np.sqrt(2.0)*np.sqrt(2.0))))))) + assert_equal(y,1) + +class TestBasicStats(TestCase): + """ W.II.C. Compute basic statistic on all the variables. + + The means should be the fifth value of all the variables (case FIVE). + The standard deviations should be "undefined" or missing for MISS, + 0 for ZERO, and 2.738612788 (times 10 to a power) for all the other variables. + II. C. Basic Statistics + """ + + dprec = np.finfo(np.float64).precision + + # Really need to write these tests to handle missing values properly + def test_tmeanX(self): + y = stats.tmean(X, (2, 8), (True, True)) + assert_approx_equal(y, 5.0, significant=TestBasicStats.dprec) + + def test_tvarX(self): + y = stats.tvar(X, (2, 8), (True, True)) + assert_approx_equal(y, 4.6666666666666661, + significant=TestBasicStats.dprec) + + def test_tstdX(self): + y = stats.tstd(X, (2, 8), (True, True)) + assert_approx_equal(y, 2.1602468994692865, + significant=TestBasicStats.dprec) + + + +class TestNanFunc(TestCase): + def __init__(self, *args, **kw): + TestCase.__init__(self, *args, **kw) + self.X = X.copy() + + self.Xall = X.copy() + self.Xall[:] = np.nan + + self.Xsome = X.copy() + self.Xsomet = X.copy() + self.Xsome[0] = np.nan + self.Xsomet = self.Xsomet[1:] + + def test_nanmean_none(self): + """Check nanmean when no values are nan.""" + m = stats.nanmean(X) + assert_approx_equal(m, X[4]) + + def test_nanmean_some(self): + """Check nanmean when some values only are nan.""" + m = stats.nanmean(self.Xsome) + assert_approx_equal(m, 5.5) + + def test_nanmean_all(self): + """Check nanmean when all values are nan.""" + olderr = np.seterr(all='ignore') + try: + m = stats.nanmean(self.Xall) + finally: + np.seterr(**olderr) + assert_(np.isnan(m)) + + def test_nanstd_none(self): + """Check nanstd when no values are nan.""" + s = stats.nanstd(self.X) + assert_approx_equal(s, np.std(self.X, ddof=1)) + + def test_nanstd_some(self): + """Check nanstd when some values only are nan.""" + s = stats.nanstd(self.Xsome) + assert_approx_equal(s, np.std(self.Xsomet, ddof=1)) + + def test_nanstd_all(self): + """Check nanstd when all values are nan.""" + olderr = np.seterr(all='ignore') + try: + s = stats.nanstd(self.Xall) + finally: + np.seterr(**olderr) + assert_(np.isnan(s)) + + def test_nanstd_negative_axis(self): + x = np.array([1, 2, 3]) + assert_equal(stats.nanstd(x, -1), 1) + + def test_nanmedian_none(self): + """Check nanmedian when no values are nan.""" + m = stats.nanmedian(self.X) + assert_approx_equal(m, np.median(self.X)) + + def test_nanmedian_some(self): + """Check nanmedian when some values only are nan.""" + m = stats.nanmedian(self.Xsome) + assert_approx_equal(m, np.median(self.Xsomet)) + + def test_nanmedian_all(self): + """Check nanmedian when all values are nan.""" + m = stats.nanmedian(self.Xall) + assert_(np.isnan(m)) + + def test_nanmedian_scalars(self): + """Check nanmedian for scalar inputs. See ticket #1098.""" + assert_equal(stats.nanmedian(1), np.median(1)) + assert_equal(stats.nanmedian(True), np.median(True)) + assert_equal(stats.nanmedian(np.array(1)), np.median(np.array(1))) + assert_equal(stats.nanmedian(np.nan), np.median(np.nan)) + + +class TestCorrPearsonr(TestCase): + """ W.II.D. Compute a correlation matrix on all the variables. + + All the correlations, except for ZERO and MISS, shoud be exactly 1. + ZERO and MISS should have undefined or missing correlations with the + other variables. The same should go for SPEARMAN corelations, if + your program has them. + """ + def test_pXX(self): + y = stats.pearsonr(X,X) + r = y[0] + assert_approx_equal(r,1.0) + + def test_pXBIG(self): + y = stats.pearsonr(X,BIG) + r = y[0] + assert_approx_equal(r,1.0) + + def test_pXLITTLE(self): + y = stats.pearsonr(X,LITTLE) + r = y[0] + assert_approx_equal(r,1.0) + + def test_pXHUGE(self): + y = stats.pearsonr(X,HUGE) + r = y[0] + assert_approx_equal(r,1.0) + + def test_pXTINY(self): + y = stats.pearsonr(X,TINY) + r = y[0] + assert_approx_equal(r,1.0) + + def test_pXROUND(self): + y = stats.pearsonr(X,ROUND) + r = y[0] + assert_approx_equal(r,1.0) + + def test_pBIGBIG(self): + y = stats.pearsonr(BIG,BIG) + r = y[0] + assert_approx_equal(r,1.0) + + def test_pBIGLITTLE(self): + y = stats.pearsonr(BIG,LITTLE) + r = y[0] + assert_approx_equal(r,1.0) + + def test_pBIGHUGE(self): + y = stats.pearsonr(BIG,HUGE) + r = y[0] + assert_approx_equal(r,1.0) + + def test_pBIGTINY(self): + y = stats.pearsonr(BIG,TINY) + r = y[0] + assert_approx_equal(r,1.0) + + def test_pBIGROUND(self): + y = stats.pearsonr(BIG,ROUND) + r = y[0] + assert_approx_equal(r,1.0) + + def test_pLITTLELITTLE(self): + y = stats.pearsonr(LITTLE,LITTLE) + r = y[0] + assert_approx_equal(r,1.0) + + def test_pLITTLEHUGE(self): + y = stats.pearsonr(LITTLE,HUGE) + r = y[0] + assert_approx_equal(r,1.0) + + def test_pLITTLETINY(self): + y = stats.pearsonr(LITTLE,TINY) + r = y[0] + assert_approx_equal(r,1.0) + + def test_pLITTLEROUND(self): + y = stats.pearsonr(LITTLE,ROUND) + r = y[0] + assert_approx_equal(r,1.0) + + def test_pHUGEHUGE(self): + y = stats.pearsonr(HUGE,HUGE) + r = y[0] + assert_approx_equal(r,1.0) + + def test_pHUGETINY(self): + y = stats.pearsonr(HUGE,TINY) + r = y[0] + assert_approx_equal(r,1.0) + + def test_pHUGEROUND(self): + y = stats.pearsonr(HUGE,ROUND) + r = y[0] + assert_approx_equal(r,1.0) + + def test_pTINYTINY(self): + y = stats.pearsonr(TINY,TINY) + r = y[0] + assert_approx_equal(r,1.0) + + def test_pTINYROUND(self): + y = stats.pearsonr(TINY,ROUND) + r = y[0] + assert_approx_equal(r,1.0) + + def test_pROUNDROUND(self): + y = stats.pearsonr(ROUND,ROUND) + r = y[0] + assert_approx_equal(r,1.0) + + def test_r_exactly_pos1(self): + a = arange(3.0) + b = a + r, prob = stats.pearsonr(a,b) + assert_equal(r, 1.0) + assert_equal(prob, 0.0) + + def test_r_exactly_neg1(self): + a = arange(3.0) + b = -a + r, prob = stats.pearsonr(a,b) + assert_equal(r, -1.0) + assert_equal(prob, 0.0) + +class TestFisherExact(TestCase): + """Some tests to show that fisher_exact() works correctly. + + Note that in SciPy 0.9.0 this was not working well for large numbers due to + inaccuracy of the hypergeom distribution (see #1218). Fixed now. + + Also note that R and Scipy have different argument formats for their + hypergeometric distribution functions. + + R: + > phyper(18999, 99000, 110000, 39000, lower.tail = FALSE) + [1] 1.701815e-09 + """ + def test_basic(self): + fisher_exact = stats.fisher_exact + + res = fisher_exact([[14500, 20000], [30000, 40000]])[1] + assert_approx_equal(res, 0.01106, significant=4) + res = fisher_exact([[100, 2], [1000, 5]])[1] + assert_approx_equal(res, 0.1301, significant=4) + res = fisher_exact([[2, 7], [8, 2]])[1] + assert_approx_equal(res, 0.0230141, significant=6) + res = fisher_exact([[5, 1], [10, 10]])[1] + assert_approx_equal(res, 0.1973244, significant=6) + res = fisher_exact([[5, 15], [20, 20]])[1] + assert_approx_equal(res, 0.0958044, significant=6) + res = fisher_exact([[5, 16], [20, 25]])[1] + assert_approx_equal(res, 0.1725862, significant=6) + res = fisher_exact([[10, 5], [10, 1]])[1] + assert_approx_equal(res, 0.1973244, significant=6) + res = fisher_exact([[5, 0], [1, 4]])[1] + assert_approx_equal(res, 0.04761904, significant=6) + res = fisher_exact([[0, 1], [3, 2]])[1] + assert_approx_equal(res, 1.0) + res = fisher_exact([[0, 2], [6, 4]])[1] + assert_approx_equal(res, 0.4545454545) + res = fisher_exact([[2, 7], [8, 2]]) + assert_approx_equal(res[1], 0.0230141, significant=6) + assert_approx_equal(res[0], 4.0 / 56) + + def test_precise(self): + fisher_exact = stats.fisher_exact + + # results from R + # + # R defines oddsratio differently (see Notes section of fisher_exact + # docstring), so those will not match. We leave them in anyway, in + # case they will be useful later on. We test only the p-value. + tablist = [ + ([[100, 2], [1000, 5]], (2.505583993422285e-001, 1.300759363430016e-001)), + ([[2, 7], [8, 2]], (8.586235135736206e-002, 2.301413756522114e-002)), + ([[5, 1], [10, 10]], (4.725646047336584e+000, 1.973244147157190e-001)), + ([[5, 15], [20, 20]], (3.394396617440852e-001, 9.580440012477637e-002)), + ([[5, 16], [20, 25]], (3.960558326183334e-001, 1.725864953812994e-001)), + ([[10, 5], [10, 1]], (2.116112781158483e-001, 1.973244147157190e-001)), + ([[10, 5], [10, 0]], (0.000000000000000e+000, 6.126482213438734e-002)), + ([[5, 0], [1, 4]], (np.inf, 4.761904761904762e-002)), + ([[0, 5], [1, 4]], (0.000000000000000e+000, 1.000000000000000e+000)), + ([[5, 1], [0, 4]], (np.inf, 4.761904761904758e-002)), + ([[0, 1], [3, 2]], (0.000000000000000e+000, 1.000000000000000e+000)) + ] + for table, res_r in tablist: + res = stats.fisher_exact(np.asarray(table)) + np.testing.assert_almost_equal(res[1], res_r[1], decimal=11, + verbose=True) + + @dec.slow + def test_large_numbers(self): + # Test with some large numbers. Regression test for #1401 + pvals = [5.56e-11, 2.666e-11, 1.363e-11] # from R + for pval, num in zip(pvals, [75, 76, 77]): + res = stats.fisher_exact([[17704, 496], [1065, num]])[1] + assert_approx_equal(res, pval, significant=4) + + res = stats.fisher_exact([[18000, 80000], [20000, 90000]])[1] + assert_approx_equal(res, 0.2751, significant=4) + + def test_raises(self): + # test we raise an error for wrong shape of input. + assert_raises(ValueError, stats.fisher_exact, + np.arange(6).reshape(2, 3)) + + def test_row_or_col_zero(self): + tables = ([[0, 0], [5, 10]], + [[5, 10], [0, 0]], + [[0, 5], [0, 10]], + [[5, 0], [10, 0]]) + for table in tables: + oddsratio, pval = stats.fisher_exact(table) + assert_equal(pval, 1.0) + assert_equal(oddsratio, np.nan) + + def test_less_greater(self): + tables = ( + # Some tables to compare with R: + [[2, 7], [8, 2]], + [[200, 7], [8, 300]], + [[28, 21], [6, 1957]], + [[190, 800], [200, 900]], + # Some tables with simple exact values + # (includes regression test for ticket #1568): + [[0, 2], [3, 0]], + [[1, 1], [2, 1]], + [[2, 0], [1, 2]], + [[0, 1], [2, 3]], + [[1, 0], [1, 4]], + ) + pvals = ( + # from R: + [0.018521725952066501, 0.9990149169715733], + [1.0, 2.0056578803889148e-122], + [1.0, 5.7284374608319831e-44], + [0.7416227, 0.2959826], + # Exact: + [0.1, 1.0], + [0.7, 0.9], + [1.0, 0.3], + [2./3, 1.0], + [1.0, 1./3], + ) + for table, pval in zip(tables, pvals): + res = [] + res.append(stats.fisher_exact(table, alternative="less")[1]) + res.append(stats.fisher_exact(table, alternative="greater")[1]) + assert_allclose(res, pval, atol=0, rtol=1e-7) + + +class TestCorrSpearmanr(TestCase): + """ W.II.D. Compute a correlation matrix on all the variables. + + All the correlations, except for ZERO and MISS, shoud be exactly 1. + ZERO and MISS should have undefined or missing correlations with the + other variables. The same should go for SPEARMAN corelations, if + your program has them. + """ + def test_sXX(self): + y = stats.spearmanr(X,X) + r = y[0] + assert_approx_equal(r,1.0) + + def test_sXBIG(self): + y = stats.spearmanr(X,BIG) + r = y[0] + assert_approx_equal(r,1.0) + + def test_sXLITTLE(self): + y = stats.spearmanr(X,LITTLE) + r = y[0] + assert_approx_equal(r,1.0) + + def test_sXHUGE(self): + y = stats.spearmanr(X,HUGE) + r = y[0] + assert_approx_equal(r,1.0) + + def test_sXTINY(self): + y = stats.spearmanr(X,TINY) + r = y[0] + assert_approx_equal(r,1.0) + + def test_sXROUND(self): + y = stats.spearmanr(X,ROUND) + r = y[0] + assert_approx_equal(r,1.0) + + def test_sBIGBIG(self): + y = stats.spearmanr(BIG,BIG) + r = y[0] + assert_approx_equal(r,1.0) + + def test_sBIGLITTLE(self): + y = stats.spearmanr(BIG,LITTLE) + r = y[0] + assert_approx_equal(r,1.0) + + def test_sBIGHUGE(self): + y = stats.spearmanr(BIG,HUGE) + r = y[0] + assert_approx_equal(r,1.0) + + def test_sBIGTINY(self): + y = stats.spearmanr(BIG,TINY) + r = y[0] + assert_approx_equal(r,1.0) + + def test_sBIGROUND(self): + y = stats.spearmanr(BIG,ROUND) + r = y[0] + assert_approx_equal(r,1.0) + + def test_sLITTLELITTLE(self): + y = stats.spearmanr(LITTLE,LITTLE) + r = y[0] + assert_approx_equal(r,1.0) + + def test_sLITTLEHUGE(self): + y = stats.spearmanr(LITTLE,HUGE) + r = y[0] + assert_approx_equal(r,1.0) + + def test_sLITTLETINY(self): + y = stats.spearmanr(LITTLE,TINY) + r = y[0] + assert_approx_equal(r,1.0) + + def test_sLITTLEROUND(self): + y = stats.spearmanr(LITTLE,ROUND) + r = y[0] + assert_approx_equal(r,1.0) + + def test_sHUGEHUGE(self): + y = stats.spearmanr(HUGE,HUGE) + r = y[0] + assert_approx_equal(r,1.0) + + def test_sHUGETINY(self): + y = stats.spearmanr(HUGE,TINY) + r = y[0] + assert_approx_equal(r,1.0) + + def test_sHUGEROUND(self): + y = stats.spearmanr(HUGE,ROUND) + r = y[0] + assert_approx_equal(r,1.0) + + def test_sTINYTINY(self): + y = stats.spearmanr(TINY,TINY) + r = y[0] + assert_approx_equal(r,1.0) + + def test_sTINYROUND(self): + y = stats.spearmanr(TINY,ROUND) + r = y[0] + assert_approx_equal(r,1.0) + + def test_sROUNDROUND(self): + y = stats.spearmanr(ROUND,ROUND) + r = y[0] + assert_approx_equal(r,1.0) + +class TestCorrSpearmanrTies(TestCase): + """Some tests of tie-handling by the spearmanr function.""" + + def test_tie1(self): + # Data + x = [1.0, 2.0, 3.0, 4.0] + y = [1.0, 2.0, 2.0, 3.0] + # Ranks of the data, with tie-handling. + xr = [1.0, 2.0, 3.0, 4.0] + yr = [1.0, 2.5, 2.5, 4.0] + # Result of spearmanr should be the same as applying + # pearsonr to the ranks. + sr = stats.spearmanr(x, y) + pr = stats.pearsonr(xr, yr) + assert_almost_equal(sr, pr) + + +## W.II.E. Tabulate X against X, using BIG as a case weight. The values +## should appear on the diagonal and the total should be 899999955. +## If the table cannot hold these values, forget about working with +## census data. You can also tabulate HUGE against TINY. There is no +## reason a tabulation program should not be able to distinguish +## different values regardless of their magnitude. + +### I need to figure out how to do this one. + + +##@dec.knownfailureif(sys.version[:3] < '2.5', "Can't index array with np.int64") +@dec.slow +def test_kendalltau(): + # with some ties + x1 = [12, 2, 1, 12, 2] + x2 = [1, 4, 7, 1, 0] + expected = (-0.47140452079103173, 0.24821309157521476) + res = stats.kendalltau(x1, x2) + assert_approx_equal(res[0], expected[0]) + assert_approx_equal(res[1], expected[1]) + + # with only ties in one or both inputs + assert_(np.all(np.isnan(stats.kendalltau([2,2,2], [2,2,2])))) + assert_(np.all(np.isnan(stats.kendalltau([2,0,2], [2,2,2])))) + assert_(np.all(np.isnan(stats.kendalltau([2,2,2], [2,0,2])))) + + # check two different sort methods + assert_approx_equal(stats.kendalltau(x1, x2, initial_lexsort=False)[1], + stats.kendalltau(x1, x2, initial_lexsort=True)[1]) + + # and with larger arrays + np.random.seed(7546) + x = np.array([np.random.normal(loc=1, scale=1, size=500), + np.random.normal(loc=1, scale=1, size=500)]) + corr = [[1.0, 0.3], + [0.3, 1.0]] + x = np.dot(np.linalg.cholesky(corr), x) + expected = (0.19291382765531062, 1.1337108207276285e-10) + res = stats.kendalltau(x[0], x[1]) + assert_approx_equal(res[0], expected[0]) + assert_approx_equal(res[1], expected[1]) + + # and do we get a tau of 1 for identical inputs? + assert_approx_equal(stats.kendalltau([1,1,2], [1,1,2])[0], 1.0) + + +class TestRegression(TestCase): + def test_linregressBIGX(self): + """ W.II.F. Regress BIG on X. + + The constant should be 99999990 and the regression coefficient should be 1. + """ + y = stats.linregress(X,BIG) + intercept = y[1] + r=y[2] + assert_almost_equal(intercept,99999990) + assert_almost_equal(r,1.0) + +## W.IV.A. Take the NASTY dataset above. Use the variable X as a +## basis for computing polynomials. Namely, compute X1=X, X2=X*X, +## X3=X*X*X, and so on up to 9 products. Use the algebraic +## transformation language within the statistical package itself. You +## will end up with 9 variables. Now regress X1 on X2-X9 (a perfect +## fit). If the package balks (singular or roundoff error messages), +## try X1 on X2-X8, and so on. Most packages cannot handle more than +## a few polynomials. +## Scipy's stats.py does not seem to handle multiple linear regression +## The datasets X1 . . X9 are at the top of the file. + + + def test_regressXX(self): + """ W.IV.B. Regress X on X. + + The constant should be exactly 0 and the regression coefficient should be 1. + This is a perfectly valid regression. The program should not complain. + """ + y = stats.linregress(X,X) + intercept = y[1] + r=y[2] + assert_almost_equal(intercept,0.0) + assert_almost_equal(r,1.0) +## W.IV.C. Regress X on BIG and LITTLE (two predictors). The program +## should tell you that this model is "singular" because BIG and +## LITTLE are linear combinations of each other. Cryptic error +## messages are unacceptable here. Singularity is the most +## fundamental regression error. +### Need to figure out how to handle multiple linear regression. Not obvious + + def test_regressZEROX(self): + """ W.IV.D. Regress ZERO on X. + + The program should inform you that ZERO has no variance or it should + go ahead and compute the regression and report a correlation and + total sum of squares of exactly 0. + """ + y = stats.linregress(X,ZERO) + intercept = y[1] + r=y[2] + assert_almost_equal(intercept,0.0) + assert_almost_equal(r,0.0) + + def test_regress_simple(self): + """Regress a line with sinusoidal noise.""" + x = np.linspace(0, 100, 100) + y = 0.2 * np.linspace(0, 100, 100) + 10 + y += np.sin(np.linspace(0, 20, 100)) + + res = stats.linregress(x, y) + assert_almost_equal(res[4], 2.3957814497838803e-3) #4.3609875083149268e-3) + + def test_regress_simple_onearg_rows(self): + """Regress a line with sinusoidal noise, with a single input of shape + (2, N). + """ + x = np.linspace(0, 100, 100) + y = 0.2 * np.linspace(0, 100, 100) + 10 + y += np.sin(np.linspace(0, 20, 100)) + rows = np.vstack((x, y)) + + res = stats.linregress(rows) + assert_almost_equal(res[4], 2.3957814497838803e-3) #4.3609875083149268e-3) + + def test_regress_simple_onearg_cols(self): + """Regress a line with sinusoidal noise, with a single input of shape + (N, 2). + """ + x = np.linspace(0, 100, 100) + y = 0.2 * np.linspace(0, 100, 100) + 10 + y += np.sin(np.linspace(0, 20, 100)) + cols = np.hstack((np.expand_dims(x, 1), np.expand_dims(y, 1))) + + res = stats.linregress(cols) + assert_almost_equal(res[4], 2.3957814497838803e-3) #4.3609875083149268e-3) + + def test_regress_shape_error(self): + """Check that a single input argument to linregress with wrong shape + results in a ValueError.""" + assert_raises(ValueError, stats.linregress, np.ones((3, 3))) + + def test_linregress(self): + '''compared with multivariate ols with pinv''' + x = np.arange(11) + y = np.arange(5,16) + y[[(1),(-2)]] -= 1 + y[[(0),(-1)]] += 1 + + res = (1.0, 5.0, 0.98229948625750, 7.45259691e-008, 0.063564172616372733) + assert_array_almost_equal(stats.linregress(x,y),res,decimal=14) + + def test_regress_simple_negative_cor(self): + # If the slope of the regression is negative the factor R tend to -1 not 1. + # Sometimes rounding errors makes it < -1 leading to stderr being NaN + a, n = 1e-71, 100000 + x = np.linspace(a, 2 * a, n) + y = np.linspace(2 * a, a, n) + stats.linregress(x, y) + res = stats.linregress(x, y) + assert_(res[2] >= -1) # propagated numerical errors were not corrected + assert_almost_equal(res[2], -1) # perfect negative correlation case + assert_(not np.isnan(res[4])) # stderr should stay finite + + +class TestHistogram(TestCase): + """ Tests that histogram works as it should, and keeps old behaviour + """ + # what is untested: + # - multidimensional arrays (since 'a' is ravel'd as the first line in the method) + # - very large arrays + # - Nans, Infs, empty and otherwise bad inputs + + # sample arrays to test the histogram with + low_values = np.array([0.2, 0.3, 0.4, 0.5, 0.5, 0.6, 0.7, 0.8, 0.9, 1.1, 1.2], + dtype=float) # 11 values + high_range = np.array([2, 3, 4, 2, 21, 32, 78, 95, 65, 66, 66, 66, 66, 4], + dtype=float) # 14 values + low_range = np.array([2, 3, 3, 2, 3, 2.4, 2.1, 3.1, 2.9, 2.6, 2.7, 2.8, 2.2, 2.001], + dtype=float) # 14 values + few_values = np.array([2.0, 3.0, -1.0, 0.0], dtype=float) # 4 values + + def test_simple(self): + """ Tests that each of the tests works as expected with default params + """ + # basic tests, with expected results (no weighting) + # results taken from the previous (slower) version of histogram + basic_tests = ((self.low_values, (np.array([ 1., 1., 1., 2., 2., + 1., 1., 0., 1., 1.]), + 0.14444444444444446, 0.11111111111111112, 0)), + (self.high_range, (np.array([ 5., 0., 1., 1., 0., + 0., 5., 1., 0., 1.]), + -3.1666666666666661, 10.333333333333332, 0)), + (self.low_range, (np.array([ 3., 1., 1., 1., 0., 1., + 1., 2., 3., 1.]), + 1.9388888888888889, 0.12222222222222223, 0)), + (self.few_values, (np.array([ 1., 0., 1., 0., 0., 0., + 0., 1., 0., 1.]), + -1.2222222222222223, 0.44444444444444448, 0)), + ) + for inputs, expected_results in basic_tests: + given_results = stats.histogram(inputs) + assert_array_almost_equal(expected_results[0], given_results[0], + decimal=2) + for i in range(1, 4): + assert_almost_equal(expected_results[i], given_results[i], + decimal=2) + + def test_weighting(self): + """ Tests that weights give expected histograms + """ + # basic tests, with expected results, given a set of weights + # weights used (first n are used for each test, where n is len of array) (14 values) + weights = np.array([1., 3., 4.5, 0.1, -1.0, 0.0, 0.3, 7.0, 103.2, 2, 40, 0, 0, 1]) + # results taken from the numpy version of histogram + basic_tests = ((self.low_values, (np.array([ 4.0, 0.0, 4.5, -0.9, 0.0, + 0.3,110.2, 0.0, 0.0, 42.0]), + 0.2, 0.1, 0)), + (self.high_range, (np.array([ 9.6, 0. , -1. , 0. , 0. , + 0. ,145.2, 0. , 0.3, 7. ]), + 2.0, 9.3, 0)), + (self.low_range, (np.array([ 2.4, 0. , 0. , 0. , 0. , + 2. , 40. , 0. , 103.2, 13.5]), + 2.0, 0.11, 0)), + (self.few_values, (np.array([ 4.5, 0. , 0.1, 0. , 0. , 0. , + 0. , 1. , 0. , 3. ]), + -1., 0.4, 0)), + + ) + for inputs, expected_results in basic_tests: + # use the first lot of weights for test + # default limits given to reproduce output of numpy's test better + given_results = stats.histogram(inputs, defaultlimits=(inputs.min(), + inputs.max()), + weights=weights[:len(inputs)]) + assert_array_almost_equal(expected_results[0], given_results[0], + decimal=2) + for i in range(1, 4): + assert_almost_equal(expected_results[i], given_results[i], + decimal=2) + + def test_reduced_bins(self): + """ Tests that reducing the number of bins produces expected results + """ + # basic tests, with expected results (no weighting), + # except number of bins is halved to 5 + # results taken from the previous (slower) version of histogram + basic_tests = ((self.low_values, (np.array([ 2., 3., 3., 1., 2.]), + 0.075000000000000011, 0.25, 0)), + (self.high_range, (np.array([ 5., 2., 0., 6., 1.]), + -9.625, 23.25, 0)), + (self.low_range, (np.array([ 4., 2., 1., 3., 4.]), + 1.8625, 0.27500000000000002, 0)), + (self.few_values, (np.array([ 1., 1., 0., 1., 1.]), + -1.5, 1.0, 0)), + ) + for inputs, expected_results in basic_tests: + given_results = stats.histogram(inputs, numbins=5) + assert_array_almost_equal(expected_results[0], given_results[0], + decimal=2) + for i in range(1, 4): + assert_almost_equal(expected_results[i], given_results[i], + decimal=2) + + def test_increased_bins(self): + """ Tests that increasing the number of bins produces expected results + """ + # basic tests, with expected results (no weighting), + # except number of bins is double to 20 + # results taken from the previous (slower) version of histogram + basic_tests = ((self.low_values, (np.array([ 1., 0., 1., 0., 1., + 0., 2., 0., 1., 0., + 1., 1., 0., 1., 0., + 0., 0., 1., 0., 1.]), + 0.1736842105263158, 0.052631578947368418, 0)), + (self.high_range, (np.array([ 5., 0., 0., 0., 1., + 0., 1., 0., 0., 0., + 0., 0., 0., 5., 0., + 0., 1., 0., 0., 1.]), + -0.44736842105263142, 4.8947368421052628, 0)), + (self.low_range, (np.array([ 3., 0., 1., 1., 0., 0., + 0., 1., 0., 0., 1., 0., + 1., 0., 1., 0., 1., 3., + 0., 1.]), + 1.9710526315789474, 0.057894736842105263, 0)), + (self.few_values, (np.array([ 1., 0., 0., 0., 0., 1., + 0., 0., 0., 0., 0., 0., + 0., 0., 1., 0., 0., 0., + 0., 1.]), + -1.1052631578947367, 0.21052631578947367, 0)), + ) + for inputs, expected_results in basic_tests: + given_results = stats.histogram(inputs, numbins=20) + assert_array_almost_equal(expected_results[0], given_results[0], + decimal=2) + for i in range(1, 4): + assert_almost_equal(expected_results[i], given_results[i], + decimal=2) + + +def test_cumfreq(): + x = [1, 4, 2, 1, 3, 1] + cumfreqs, lowlim, binsize, extrapoints = stats.cumfreq(x, numbins=4) + assert_array_almost_equal(cumfreqs, np.array([ 3., 4., 5., 6.])) + cumfreqs, lowlim, binsize, extrapoints = stats.cumfreq(x, numbins=4, + defaultreallimits=(1.5, 5)) + assert_(extrapoints==3) + + +def test_relfreq(): + a = np.array([1, 4, 2, 1, 3, 1]) + relfreqs, lowlim, binsize, extrapoints = stats.relfreq(a, numbins=4) + assert_array_almost_equal(relfreqs, array([0.5, 0.16666667, 0.16666667, 0.16666667])) + + # check array_like input is accepted + relfreqs2, lowlim, binsize, extrapoints = stats.relfreq([1, 4, 2, 1, 3, 1], numbins=4) + assert_array_almost_equal(relfreqs, relfreqs2) + + +# Utility + +def compare_results(res,desired): + for i in range(len(desired)): + assert_array_equal(res[i],desired[i]) + + +################################################## +### Test for sum + +class TestGMean(TestCase): + + def test_1D_list(self): + a = (1,2,3,4) + actual= stats.gmean(a) + desired = power(1*2*3*4,1./4.) + assert_almost_equal(actual, desired,decimal=14) + + desired1 = stats.gmean(a,axis=-1) + assert_almost_equal(actual, desired1, decimal=14) + + def test_1D_array(self): + a = array((1,2,3,4), float32) + actual= stats.gmean(a) + desired = power(1*2*3*4,1./4.) + assert_almost_equal(actual, desired, decimal=7) + + desired1 = stats.gmean(a,axis=-1) + assert_almost_equal(actual, desired1, decimal=7) + + def test_2D_array_default(self): + a = array(((1,2,3,4), + (1,2,3,4), + (1,2,3,4))) + actual= stats.gmean(a) + desired = array((1,2,3,4)) + assert_array_almost_equal(actual, desired, decimal=14) + + desired1 = stats.gmean(a,axis=0) + assert_array_almost_equal(actual, desired1, decimal=14) + + def test_2D_array_dim1(self): + a = array(((1,2,3,4), + (1,2,3,4), + (1,2,3,4))) + actual= stats.gmean(a, axis=1) + v = power(1*2*3*4,1./4.) + desired = array((v,v,v)) + assert_array_almost_equal(actual, desired, decimal=14) + + def test_large_values(self): + a = array([1e100, 1e200, 1e300]) + actual = stats.gmean(a) + assert_approx_equal(actual, 1e200, significant=14) + +class TestHMean(TestCase): + def test_1D_list(self): + a = (1,2,3,4) + actual= stats.hmean(a) + desired = 4. / (1./1 + 1./2 + 1./3 + 1./4) + assert_almost_equal(actual, desired, decimal=14) + + desired1 = stats.hmean(array(a),axis=-1) + assert_almost_equal(actual, desired1, decimal=14) + def test_1D_array(self): + a = array((1,2,3,4), float64) + actual= stats.hmean(a) + desired = 4. / (1./1 + 1./2 + 1./3 + 1./4) + assert_almost_equal(actual, desired, decimal=14) + + desired1 = stats.hmean(a,axis=-1) + assert_almost_equal(actual, desired1, decimal=14) + + def test_2D_array_default(self): + a = array(((1,2,3,4), + (1,2,3,4), + (1,2,3,4))) + actual = stats.hmean(a) + desired = array((1.,2.,3.,4.)) + assert_array_almost_equal(actual, desired, decimal=14) + + actual1 = stats.hmean(a,axis=0) + assert_array_almost_equal(actual1, desired, decimal=14) + + def test_2D_array_dim1(self): + a = array(((1,2,3,4), + (1,2,3,4), + (1,2,3,4))) + + v = 4. / (1./1 + 1./2 + 1./3 + 1./4) + desired1 = array((v,v,v)) + actual1 = stats.hmean(a, axis=1) + assert_array_almost_equal(actual1, desired1, decimal=14) + +@dec.slow +class TestScoreatpercentile(TestCase): + def setUp(self): + self.a1 = [3, 4, 5, 10, -3, -5, 6] + self.a2 = [3, -6, -2, 8, 7, 4, 2, 1] + self.a3 = [3., 4, 5, 10, -3, -5, -6, 7.0] + + def test_basic(self): + x = arange(8) * 0.5 + assert_equal(stats.scoreatpercentile(x, 0), 0.) + assert_equal(stats.scoreatpercentile(x, 100), 3.5) + assert_equal(stats.scoreatpercentile(x, 50), 1.75) + + def test_2D(self): + x = array([[1, 1, 1], + [1, 1, 1], + [4, 4, 3], + [1, 1, 1], + [1, 1, 1]]) + assert_array_equal(stats.scoreatpercentile(x, 50), [1, 1, 1]) + + def test_fraction(self): + scoreatperc = stats.scoreatpercentile + + # Test defaults + assert_equal(scoreatperc(list(range(10)), 50), 4.5) + assert_equal(scoreatperc(list(range(10)), 50, (2,7)), 4.5) + assert_equal(scoreatperc(list(range(100)), 50, limit=(1, 8)), 4.5) + assert_equal(scoreatperc(np.array([1, 10 ,100]), 50, (10,100)), 55) + assert_equal(scoreatperc(np.array([1, 10 ,100]), 50, (1,10)), 5.5) + + # explicitly specify interpolation_method 'fraction' (the default) + assert_equal(scoreatperc(list(range(10)), 50, interpolation_method='fraction'), + 4.5) + assert_equal(scoreatperc(list(range(10)), 50, limit=(2, 7), + interpolation_method='fraction'), + 4.5) + assert_equal(scoreatperc(list(range(100)), 50, limit=(1, 8), + interpolation_method='fraction'), + 4.5) + assert_equal(scoreatperc(np.array([1, 10 ,100]), 50, (10, 100), + interpolation_method='fraction'), + 55) + assert_equal(scoreatperc(np.array([1, 10 ,100]), 50, (1,10), + interpolation_method='fraction'), + 5.5) + + def test_lower_higher(self): + scoreatperc = stats.scoreatpercentile + + # interpolation_method 'lower'/'higher' + assert_equal(scoreatperc(list(range(10)), 50, + interpolation_method='lower'), 4) + assert_equal(scoreatperc(list(range(10)), 50, + interpolation_method='higher'), 5) + assert_equal(scoreatperc(list(range(10)), 50, (2,7), + interpolation_method='lower'), 4) + assert_equal(scoreatperc(list(range(10)), 50, limit=(2,7), + interpolation_method='higher'), 5) + assert_equal(scoreatperc(list(range(100)), 50, (1,8), + interpolation_method='lower'), 4) + assert_equal(scoreatperc(list(range(100)), 50, (1,8), + interpolation_method='higher'), 5) + assert_equal(scoreatperc(np.array([1, 10, 100]), 50, (10, 100), + interpolation_method='lower'), 10) + assert_equal(scoreatperc(np.array([1, 10, 100]), 50, limit=(10, 100), + interpolation_method='higher'), 100) + assert_equal(scoreatperc(np.array([1, 10, 100]), 50, (1, 10), + interpolation_method='lower'), 1) + assert_equal(scoreatperc(np.array([1, 10, 100]), 50, limit=(1, 10), + interpolation_method='higher'), 10) + + def test_sequence(self): + x = arange(8) * 0.5 + assert_equal(stats.scoreatpercentile(x, [0, 100, 50]), [0, 3.5, 1.75]) + + def test_axis(self): + scoreatperc = stats.scoreatpercentile + x = arange(12).reshape(3, 4) + + assert_equal(scoreatperc(x, (25, 50, 100)), [2.75, 5.5, 11.0]) + + r0 = [[2, 3, 4, 5], [4, 5, 6, 7], [8, 9, 10, 11]] + assert_equal(scoreatperc(x, (25, 50, 100), axis=0), r0) + + r1 = [[0.75, 4.75, 8.75], [1.5, 5.5, 9.5], [3, 7, 11]] + assert_equal(scoreatperc(x, (25, 50, 100), axis=1), r1) + + def test_exception(self): + assert_raises(ValueError, stats.scoreatpercentile, [1, 2], 56, + interpolation_method='foobar') + assert_raises(ValueError, stats.scoreatpercentile, [1], 101) + assert_raises(ValueError, stats.scoreatpercentile, [1], -1) + + +class TestCMedian(TestCase): + def test_basic(self): + data = [1,2,3,1,5,3,6,4,3,2,4,3,5,2.0] + assert_almost_equal(stats.cmedian(data,5),3.2916666666666665) + assert_almost_equal(stats.cmedian(data,3),3.083333333333333) + assert_almost_equal(stats.cmedian(data),3.0020020020020022) + + +class TestMode(TestCase): + def test_basic(self): + data1 = [3,5,1,10,23,3,2,6,8,6,10,6] + vals = stats.mode(data1) + assert_almost_equal(vals[0][0],6) + assert_almost_equal(vals[1][0],3) + + +class TestVariability(TestCase): + + testcase = [1,2,3,4] + + def test_signaltonoise(self): + """ + this is not in R, so used + mean(testcase,axis=0)/(sqrt(var(testcase)*3/4)) """ + #y = stats.signaltonoise(self.shoes[0]) + #assert_approx_equal(y,4.5709967) + y = stats.signaltonoise(self.testcase) + assert_approx_equal(y,2.236067977) + + def test_sem(self): + """ + this is not in R, so used + sqrt(var(testcase)*3/4)/sqrt(3) + """ + #y = stats.sem(self.shoes[0]) + #assert_approx_equal(y,0.775177399) + y = stats.sem(self.testcase) + assert_approx_equal(y,0.6454972244) + + def test_zmap(self): + """ + not in R, so tested by using + (testcase[i]-mean(testcase,axis=0))/sqrt(var(testcase)*3/4) + """ + y = stats.zmap(self.testcase,self.testcase) + desired = ([-1.3416407864999, -0.44721359549996 , 0.44721359549996 , 1.3416407864999]) + assert_array_almost_equal(desired,y,decimal=12) + + def test_zmap_axis(self): + """Test use of 'axis' keyword in zmap.""" + x = np.array([[0.0, 0.0, 1.0, 1.0], + [1.0, 1.0, 1.0, 2.0], + [2.0, 0.0, 2.0, 0.0]]) + + t1 = 1.0/np.sqrt(2.0/3) + t2 = np.sqrt(3.)/3 + t3 = np.sqrt(2.) + + z0 = stats.zmap(x, x, axis=0) + z1 = stats.zmap(x, x, axis=1) + + z0_expected = [[-t1, -t3/2, -t3/2, 0.0], + [0.0, t3, -t3/2, t1], + [t1, -t3/2, t3, -t1]] + z1_expected = [[-1.0, -1.0, 1.0, 1.0], + [-t2, -t2, -t2, np.sqrt(3.)], + [1.0, -1.0, 1.0, -1.0]] + + assert_array_almost_equal(z0, z0_expected) + assert_array_almost_equal(z1, z1_expected) + + def test_zmap_ddof(self): + """Test use of 'ddof' keyword in zmap.""" + x = np.array([[0.0, 0.0, 1.0, 1.0], + [0.0, 1.0, 2.0, 3.0]]) + + t1 = 1.0/np.sqrt(2.0/3) + t2 = np.sqrt(3.)/3 + t3 = np.sqrt(2.) + + z = stats.zmap(x, x, axis=1, ddof=1) + + z0_expected = np.array([-0.5, -0.5, 0.5, 0.5])/(1.0/np.sqrt(3)) + z1_expected = np.array([-1.5, -0.5, 0.5, 1.5])/(np.sqrt(5./3)) + assert_array_almost_equal(z[0], z0_expected) + assert_array_almost_equal(z[1], z1_expected) + + def test_zscore(self): + """ + not in R, so tested by using + (testcase[i]-mean(testcase,axis=0))/sqrt(var(testcase)*3/4) + """ + y = stats.zscore(self.testcase) + desired = ([-1.3416407864999, -0.44721359549996 , 0.44721359549996 , 1.3416407864999]) + assert_array_almost_equal(desired,y,decimal=12) + + def test_zscore_axis(self): + """Test use of 'axis' keyword in zscore.""" + x = np.array([[0.0, 0.0, 1.0, 1.0], + [1.0, 1.0, 1.0, 2.0], + [2.0, 0.0, 2.0, 0.0]]) + + t1 = 1.0/np.sqrt(2.0/3) + t2 = np.sqrt(3.)/3 + t3 = np.sqrt(2.) + + z0 = stats.zscore(x, axis=0) + z1 = stats.zscore(x, axis=1) + + z0_expected = [[-t1, -t3/2, -t3/2, 0.0], + [0.0, t3, -t3/2, t1], + [t1, -t3/2, t3, -t1]] + z1_expected = [[-1.0, -1.0, 1.0, 1.0], + [-t2, -t2, -t2, np.sqrt(3.)], + [1.0, -1.0, 1.0, -1.0]] + + assert_array_almost_equal(z0, z0_expected) + assert_array_almost_equal(z1, z1_expected) + + def test_zscore_ddof(self): + """Test use of 'ddof' keyword in zscore.""" + x = np.array([[0.0, 0.0, 1.0, 1.0], + [0.0, 1.0, 2.0, 3.0]]) + + t1 = 1.0/np.sqrt(2.0/3) + t2 = np.sqrt(3.)/3 + t3 = np.sqrt(2.) + + z = stats.zscore(x, axis=1, ddof=1) + + z0_expected = np.array([-0.5, -0.5, 0.5, 0.5])/(1.0/np.sqrt(3)) + z1_expected = np.array([-1.5, -0.5, 0.5, 1.5])/(np.sqrt(5./3)) + assert_array_almost_equal(z[0], z0_expected) + assert_array_almost_equal(z[1], z1_expected) + + +class TestMoments(TestCase): + """ + Comparison numbers are found using R v.1.5.1 + note that length(testcase) = 4 + testmathworks comes from documentation for the + Statistics Toolbox for Matlab and can be found at both + http://www.mathworks.com/access/helpdesk/help/toolbox/stats/kurtosis.shtml + http://www.mathworks.com/access/helpdesk/help/toolbox/stats/skewness.shtml + Note that both test cases came from here. + """ + testcase = [1,2,3,4] + testmathworks = [1.165 , 0.6268, 0.0751, 0.3516, -0.6965] + def test_moment(self): + """ + mean((testcase-mean(testcase))**power,axis=0),axis=0))**power))""" + y = stats.moment(self.testcase,1) + assert_approx_equal(y,0.0,10) + y = stats.moment(self.testcase,2) + assert_approx_equal(y,1.25) + y = stats.moment(self.testcase,3) + assert_approx_equal(y,0.0) + y = stats.moment(self.testcase,4) + assert_approx_equal(y,2.5625) + def test_variation(self): + """ + variation = samplestd/mean """ +## y = stats.variation(self.shoes[0]) +## assert_approx_equal(y,21.8770668) + y = stats.variation(self.testcase) + assert_approx_equal(y,0.44721359549996, 10) + + def test_skewness(self): + """ + sum((testmathworks-mean(testmathworks,axis=0))**3,axis=0)/ + ((sqrt(var(testmathworks)*4/5))**3)/5 + """ + y = stats.skew(self.testmathworks) + assert_approx_equal(y,-0.29322304336607,10) + y = stats.skew(self.testmathworks,bias=0) + assert_approx_equal(y,-0.437111105023940,10) + y = stats.skew(self.testcase) + assert_approx_equal(y,0.0,10) + + def test_skewness_scalar(self): + """ + `skew` must return a scalar for 1-dim input + """ + assert_equal(stats.skew(arange(10)), 0.0) + + def test_kurtosis(self): + """ + sum((testcase-mean(testcase,axis=0))**4,axis=0)/((sqrt(var(testcase)*3/4))**4)/4 + sum((test2-mean(testmathworks,axis=0))**4,axis=0)/((sqrt(var(testmathworks)*4/5))**4)/5 + Set flags for axis = 0 and + fisher=0 (Pearson's defn of kurtosis for compatiability with Matlab) + """ + y = stats.kurtosis(self.testmathworks,0,fisher=0,bias=1) + assert_approx_equal(y, 2.1658856802973,10) + + # Note that MATLAB has confusing docs for the following case + # kurtosis(x,0) gives an unbiased estimate of Pearson's skewness + # kurtosis(x) gives a biased estimate of Fisher's skewness (Pearson-3) + # The MATLAB docs imply that both should give Fisher's + y = stats.kurtosis(self.testmathworks,fisher=0,bias=0) + assert_approx_equal(y, 3.663542721189047,10) + y = stats.kurtosis(self.testcase,0,0) + assert_approx_equal(y,1.64) + + def test_kurtosis_array_scalar(self): + assert_equal(type(stats.kurtosis([1,2,3])), float) + +class TestThreshold(TestCase): + def test_basic(self): + a = [-1,2,3,4,5,-1,-2] + assert_array_equal(stats.threshold(a),a) + assert_array_equal(stats.threshold(a,3,None,0), + [0,0,3,4,5,0,0]) + assert_array_equal(stats.threshold(a,None,3,0), + [-1,2,3,0,0,-1,-2]) + assert_array_equal(stats.threshold(a,2,4,0), + [0,2,3,4,0,0,0]) + +# Hypothesis test tests +class TestStudentTest(TestCase): + X1 = np.array([-1, 0, 1]) + X2 = np.array([0, 1, 2]) + T1_0 = 0 + P1_0 = 1 + T1_1 = -1.732051 + P1_1 = 0.2254033 + T1_2 = -3.464102 + P1_2 = 0.0741799 + T2_0 = 1.732051 + P2_0 = 0.2254033 + def test_onesample(self): + t, p = stats.ttest_1samp(self.X1, 0) + + assert_array_almost_equal(t, self.T1_0) + assert_array_almost_equal(p, self.P1_0) + + t, p = stats.ttest_1samp(self.X2, 0) + + assert_array_almost_equal(t, self.T2_0) + assert_array_almost_equal(p, self.P2_0) + + t, p = stats.ttest_1samp(self.X1, 1) + + assert_array_almost_equal(t, self.T1_1) + assert_array_almost_equal(p, self.P1_1) + + t, p = stats.ttest_1samp(self.X1, 2) + + assert_array_almost_equal(t, self.T1_2) + assert_array_almost_equal(p, self.P1_2) + + +def test_percentileofscore(): + pcos = stats.percentileofscore + + assert_equal(pcos([1,2,3,4,5,6,7,8,9,10],4), 40.0) + + for (kind, result) in [('mean', 35.0), + ('strict', 30.0), + ('weak', 40.0)]: + yield assert_equal, pcos(np.arange(10) + 1, + 4, kind=kind), \ + result + + # multiple - 2 + for (kind, result) in [('rank', 45.0), + ('strict', 30.0), + ('weak', 50.0), + ('mean', 40.0)]: + yield assert_equal, pcos([1,2,3,4,4,5,6,7,8,9], + 4, kind=kind), \ + result + + # multiple - 3 + assert_equal(pcos([1,2,3,4,4,4,5,6,7,8], 4), 50.0) + for (kind, result) in [('rank', 50.0), + ('mean', 45.0), + ('strict', 30.0), + ('weak', 60.0)]: + + yield assert_equal, pcos([1,2,3,4,4,4,5,6,7,8], + 4, kind=kind), \ + result + + # missing + for kind in ('rank', 'mean', 'strict', 'weak'): + yield assert_equal, pcos([1,2,3,5,6,7,8,9,10,11], + 4, kind=kind), \ + 30 + + #larger numbers + for (kind, result) in [('mean', 35.0), + ('strict', 30.0), + ('weak', 40.0)]: + yield assert_equal, \ + pcos([10, 20, 30, 40, 50, 60, 70, 80, 90, 100], 40, + kind=kind), result + + for (kind, result) in [('mean', 45.0), + ('strict', 30.0), + ('weak', 60.0)]: + yield assert_equal, \ + pcos([10, 20, 30, 40, 40, 40, 50, 60, 70, 80], + 40, kind=kind), result + + + for kind in ('rank', 'mean', 'strict', 'weak'): + yield assert_equal, \ + pcos([10, 20, 30, 50, 60, 70, 80, 90, 100, 110], + 40, kind=kind), 30.0 + + #boundaries + for (kind, result) in [('rank', 10.0), + ('mean', 5.0), + ('strict', 0.0), + ('weak', 10.0)]: + yield assert_equal, \ + pcos([10, 20, 30, 50, 60, 70, 80, 90, 100, 110], + 10, kind=kind), result + + for (kind, result) in [('rank', 100.0), + ('mean', 95.0), + ('strict', 90.0), + ('weak', 100.0)]: + yield assert_equal, \ + pcos([10, 20, 30, 50, 60, 70, 80, 90, 100, 110], + 110, kind=kind), result + + #out of bounds + for (kind, score, result) in [('rank', 200, 100.0), + ('mean', 200, 100.0), + ('mean', 0, 0.0)]: + yield assert_equal, \ + pcos([10, 20, 30, 50, 60, 70, 80, 90, 100, 110], + score, kind=kind), result + + +def test_friedmanchisquare(): + # see ticket:113 + # verified with matlab and R + #From Demsar "Statistical Comparisons of Classifiers over Multiple Data Sets" + #2006, Xf=9.28 (no tie handling, tie corrected Xf >=9.28) + x1 = [array([0.763, 0.599, 0.954, 0.628, 0.882, 0.936, 0.661, 0.583, + 0.775, 1.0, 0.94, 0.619, 0.972, 0.957]), + array([0.768, 0.591, 0.971, 0.661, 0.888, 0.931, 0.668, 0.583, + 0.838, 1.0, 0.962, 0.666, 0.981, 0.978]), + array([0.771, 0.590, 0.968, 0.654, 0.886, 0.916, 0.609, 0.563, + 0.866, 1.0, 0.965, 0.614, 0.9751, 0.946]), + array([0.798, 0.569, 0.967, 0.657, 0.898, 0.931, 0.685, 0.625, + 0.875, 1.0, 0.962, 0.669, 0.975, 0.970])] + + #From "Bioestadistica para las ciencias de la salud" Xf=18.95 p<0.001: + x2 = [array([4,3,5,3,5,3,2,5,4,4,4,3]), + array([2,2,1,2,3,1,2,3,2,1,1,3]), + array([2,4,3,3,4,3,3,4,4,1,2,1]), + array([3,5,4,3,4,4,3,3,3,4,4,4])] + + #From Jerrorl H. Zar, "Biostatistical Analysis"(example 12.6), Xf=10.68, 0.005 < p < 0.01: + #Probability from this example is inexact using Chisquare aproximation of Friedman Chisquare. + x3 = [array([7.0,9.9,8.5,5.1,10.3]), + array([5.3,5.7,4.7,3.5,7.7]), + array([4.9,7.6,5.5,2.8,8.4]), + array([8.8,8.9,8.1,3.3,9.1])] + + + assert_array_almost_equal(stats.friedmanchisquare(x1[0],x1[1],x1[2],x1[3]),(10.2283464566929, 0.0167215803284414)) + assert_array_almost_equal(stats.friedmanchisquare(x2[0],x2[1],x2[2],x2[3]),(18.9428571428571, 0.000280938375189499)) + assert_array_almost_equal(stats.friedmanchisquare(x3[0],x3[1],x3[2],x3[3]),(10.68, 0.0135882729582176)) + np.testing.assert_raises(ValueError, stats.friedmanchisquare,x3[0],x3[1]) + + # test using mstats + assert_array_almost_equal(stats.mstats.friedmanchisquare(x1[0],x1[1],x1[2],x1[3]),(10.2283464566929, 0.0167215803284414)) + # the following fails + #assert_array_almost_equal(stats.mstats.friedmanchisquare(x2[0],x2[1],x2[2],x2[3]),(18.9428571428571, 0.000280938375189499)) + assert_array_almost_equal(stats.mstats.friedmanchisquare(x3[0],x3[1],x3[2],x3[3]),(10.68, 0.0135882729582176)) + np.testing.assert_raises(ValueError,stats.mstats.friedmanchisquare,x3[0],x3[1]) + +def test_kstest(): + #from numpy.testing import assert_almost_equal + + # comparing with values from R + x = np.linspace(-1,1,9) + D,p = stats.kstest(x,'norm') + assert_almost_equal( D, 0.15865525393145705, 12) + assert_almost_equal( p, 0.95164069201518386, 1) + + x = np.linspace(-15,15,9) + D,p = stats.kstest(x,'norm') + assert_almost_equal( D, 0.44435602715924361, 15) + assert_almost_equal( p, 0.038850140086788665, 8) + + # the following tests rely on deterministicaly replicated rvs + np.random.seed(987654321) + x = stats.norm.rvs(loc=0.2, size=100) + D,p = stats.kstest(x, 'norm', mode='asymp') + assert_almost_equal( D, 0.12464329735846891, 15) + assert_almost_equal( p, 0.089444888711820769, 15) + assert_almost_equal( np.array(stats.kstest(x, 'norm', mode='asymp')), + np.array((0.12464329735846891, 0.089444888711820769)), 15) + assert_almost_equal( np.array(stats.kstest(x,'norm', alternative = 'less')), + np.array((0.12464329735846891, 0.040989164077641749)), 15) + # this 'greater' test fails with precision of decimal=14 + assert_almost_equal( np.array(stats.kstest(x,'norm', alternative = 'greater')), + np.array((0.0072115233216310994, 0.98531158590396228)), 12) + + #missing: no test that uses *args + + +def test_ks_2samp(): + #exact small sample solution + data1 = np.array([1.0,2.0]) + data2 = np.array([1.0,2.0,3.0]) + assert_almost_equal(np.array(stats.ks_2samp(data1+0.01,data2)), + np.array((0.33333333333333337, 0.99062316386915694))) + assert_almost_equal(np.array(stats.ks_2samp(data1-0.01,data2)), + np.array((0.66666666666666674, 0.42490954988801982))) + #these can also be verified graphically + assert_almost_equal( + np.array(stats.ks_2samp(np.linspace(1,100,100), + np.linspace(1,100,100)+2+0.1)), + np.array((0.030000000000000027, 0.99999999996005062))) + assert_almost_equal( + np.array(stats.ks_2samp(np.linspace(1,100,100), + np.linspace(1,100,100)+2-0.1)), + np.array((0.020000000000000018, 0.99999999999999933))) + #these are just regression tests + assert_almost_equal( + np.array(stats.ks_2samp(np.linspace(1,100,100), + np.linspace(1,100,110)+20.1)), + np.array((0.21090909090909091, 0.015880386730710221))) + assert_almost_equal( + np.array(stats.ks_2samp(np.linspace(1,100,100), + np.linspace(1,100,110)+20-0.1)), + np.array((0.20818181818181825, 0.017981441789762638))) + +def test_ttest_rel(): + #regression test + tr,pr = 0.81248591389165692, 0.41846234511362157 + tpr = ([tr,-tr],[pr,pr]) + + rvs1 = np.linspace(1,100,100) + rvs2 = np.linspace(1.01,99.989,100) + rvs1_2D = np.array([np.linspace(1,100,100), np.linspace(1.01,99.989,100)]) + rvs2_2D = np.array([np.linspace(1.01,99.989,100), np.linspace(1,100,100)]) + + t,p = stats.ttest_rel(rvs1, rvs2, axis=0) + assert_array_almost_equal([t,p],(tr,pr)) + t,p = stats.ttest_rel(rvs1_2D.T, rvs2_2D.T, axis=0) + assert_array_almost_equal([t,p],tpr) + t,p = stats.ttest_rel(rvs1_2D, rvs2_2D, axis=1) + assert_array_almost_equal([t,p],tpr) + + #test on 3 dimensions + rvs1_3D = np.dstack([rvs1_2D,rvs1_2D,rvs1_2D]) + rvs2_3D = np.dstack([rvs2_2D,rvs2_2D,rvs2_2D]) + t,p = stats.ttest_rel(rvs1_3D, rvs2_3D, axis=1) + assert_array_almost_equal(np.abs(t), tr) + assert_array_almost_equal(np.abs(p), pr) + assert_equal(t.shape, (2, 3)) + + t,p = stats.ttest_rel(np.rollaxis(rvs1_3D,2), np.rollaxis(rvs2_3D,2), axis=2) + assert_array_almost_equal(np.abs(t), tr) + assert_array_almost_equal(np.abs(p), pr) + assert_equal(t.shape, (3, 2)) + + olderr = np.seterr(all='ignore') + try: + #test zero division problem + t,p = stats.ttest_rel([0,0,0],[1,1,1]) + assert_equal((np.abs(t),p), (np.inf, 0)) + assert_equal(stats.ttest_rel([0,0,0], [0,0,0]), (np.nan, np.nan)) + + #check that nan in input array result in nan output + anan = np.array([[1,np.nan],[-1,1]]) + assert_equal(stats.ttest_ind(anan, np.zeros((2,2))),([0, np.nan], [1,np.nan])) + finally: + np.seterr(**olderr) + + +def test_ttest_ind(): + #regression test + tr = 1.0912746897927283 + pr = 0.27647818616351882 + tpr = ([tr,-tr],[pr,pr]) + + rvs2 = np.linspace(1,100,100) + rvs1 = np.linspace(5,105,100) + rvs1_2D = np.array([rvs1, rvs2]) + rvs2_2D = np.array([rvs2, rvs1]) + + t,p = stats.ttest_ind(rvs1, rvs2, axis=0) + assert_array_almost_equal([t,p],(tr,pr)) + t,p = stats.ttest_ind(rvs1_2D.T, rvs2_2D.T, axis=0) + assert_array_almost_equal([t,p],tpr) + t,p = stats.ttest_ind(rvs1_2D, rvs2_2D, axis=1) + assert_array_almost_equal([t,p],tpr) + + #test on 3 dimensions + rvs1_3D = np.dstack([rvs1_2D,rvs1_2D,rvs1_2D]) + rvs2_3D = np.dstack([rvs2_2D,rvs2_2D,rvs2_2D]) + t,p = stats.ttest_ind(rvs1_3D, rvs2_3D, axis=1) + assert_almost_equal(np.abs(t), np.abs(tr)) + assert_array_almost_equal(np.abs(p), pr) + assert_equal(t.shape, (2, 3)) + + t,p = stats.ttest_ind(np.rollaxis(rvs1_3D,2), np.rollaxis(rvs2_3D,2), axis=2) + assert_array_almost_equal(np.abs(t), np.abs(tr)) + assert_array_almost_equal(np.abs(p), pr) + assert_equal(t.shape, (3, 2)) + + olderr = np.seterr(all='ignore') + try: + #test zero division problem + t,p = stats.ttest_ind([0,0,0],[1,1,1]) + assert_equal((np.abs(t),p), (np.inf, 0)) + assert_equal(stats.ttest_ind([0,0,0], [0,0,0]), (np.nan, np.nan)) + + #check that nan in input array result in nan output + anan = np.array([[1,np.nan],[-1,1]]) + assert_equal(stats.ttest_ind(anan, np.zeros((2,2))),([0, np.nan], [1,np.nan])) + finally: + np.seterr(**olderr) + +def test_ttest_ind_with_uneq_var(): + + # check vs. R + a = (1, 2, 3) + b = (1.1, 2.9, 4.2) + pr = 0.53619490753126731 + tr = -0.68649512735572582 + t, p = stats.ttest_ind(a, b, equal_var = False) + assert_array_almost_equal([t,p], [tr, pr]) + + a = (1, 2, 3, 4) + pr = 0.84354139131608286 + tr = -0.2108663315950719 + t, p = stats.ttest_ind(a, b, equal_var = False) + assert_array_almost_equal([t,p], [tr, pr]) + + #regression test + tr = 1.0912746897927283 + tr_uneq_n = 0.66745638708050492 + pr = 0.27647831993021388 + pr_uneq_n = 0.50873585065616544 + tpr = ([tr,-tr],[pr,pr]) + + rvs3 = np.linspace(1,100, 25) + rvs2 = np.linspace(1,100,100) + rvs1 = np.linspace(5,105,100) + rvs1_2D = np.array([rvs1, rvs2]) + rvs2_2D = np.array([rvs2, rvs1]) + + t,p = stats.ttest_ind(rvs1, rvs2, axis=0, equal_var = False) + assert_array_almost_equal([t,p],(tr,pr)) + t,p = stats.ttest_ind(rvs1, rvs3, axis =0, equal_var = False) + assert_array_almost_equal([t,p], (tr_uneq_n, pr_uneq_n)) + t,p = stats.ttest_ind(rvs1_2D.T, rvs2_2D.T, axis=0, equal_var = False) + assert_array_almost_equal([t,p],tpr) + t,p = stats.ttest_ind(rvs1_2D, rvs2_2D, axis=1, equal_var = False) + assert_array_almost_equal([t,p],tpr) + + #test on 3 dimensions + rvs1_3D = np.dstack([rvs1_2D,rvs1_2D,rvs1_2D]) + rvs2_3D = np.dstack([rvs2_2D,rvs2_2D,rvs2_2D]) + t,p = stats.ttest_ind(rvs1_3D, rvs2_3D, axis=1, equal_var = False) + assert_almost_equal(np.abs(t), np.abs(tr)) + assert_array_almost_equal(np.abs(p), pr) + assert_equal(t.shape, (2, 3)) + + t,p = stats.ttest_ind(np.rollaxis(rvs1_3D,2), np.rollaxis(rvs2_3D,2), + axis=2, equal_var = False) + assert_array_almost_equal(np.abs(t), np.abs(tr)) + assert_array_almost_equal(np.abs(p), pr) + assert_equal(t.shape, (3, 2)) + + olderr = np.seterr(all='ignore') + try: + #test zero division problem + t,p = stats.ttest_ind([0,0,0],[1,1,1], equal_var = False) + assert_equal((np.abs(t),p), (np.inf, 0)) + assert_equal(stats.ttest_ind([0,0,0], [0,0,0], equal_var = False), (np.nan, np.nan)) + + #check that nan in input array result in nan output + anan = np.array([[1,np.nan],[-1,1]]) + assert_equal(stats.ttest_ind(anan, np.zeros((2,2)), equal_var = False), + ([0, np.nan], [1,np.nan])) + finally: + np.seterr(**olderr) + +def test_ttest_1samp_new(): + n1, n2, n3 = (10,15,20) + rvn1 = stats.norm.rvs(loc=5,scale=10,size=(n1,n2,n3)) + + #check multidimensional array and correct axis handling + #deterministic rvn1 and rvn2 would be better as in test_ttest_rel + t1,p1 = stats.ttest_1samp(rvn1[:,:,:], np.ones((n2,n3)),axis=0) + t2,p2 = stats.ttest_1samp(rvn1[:,:,:], 1,axis=0) + t3,p3 = stats.ttest_1samp(rvn1[:,0,0], 1) + assert_array_almost_equal(t1,t2, decimal=14) + assert_almost_equal(t1[0,0],t3, decimal=14) + assert_equal(t1.shape, (n2,n3)) + + t1,p1 = stats.ttest_1samp(rvn1[:,:,:], np.ones((n1,n3)),axis=1) + t2,p2 = stats.ttest_1samp(rvn1[:,:,:], 1,axis=1) + t3,p3 = stats.ttest_1samp(rvn1[0,:,0], 1) + assert_array_almost_equal(t1,t2, decimal=14) + assert_almost_equal(t1[0,0],t3, decimal=14) + assert_equal(t1.shape, (n1,n3)) + + t1,p1 = stats.ttest_1samp(rvn1[:,:,:], np.ones((n1,n2)),axis=2) + t2,p2 = stats.ttest_1samp(rvn1[:,:,:], 1,axis=2) + t3,p3 = stats.ttest_1samp(rvn1[0,0,:], 1) + assert_array_almost_equal(t1,t2, decimal=14) + assert_almost_equal(t1[0,0],t3, decimal=14) + assert_equal(t1.shape, (n1,n2)) + + olderr = np.seterr(all='ignore') + try: + #test zero division problem + t,p = stats.ttest_1samp([0,0,0], 1) + assert_equal((np.abs(t),p), (np.inf, 0)) + assert_equal(stats.ttest_1samp([0,0,0], 0), (np.nan, np.nan)) + + #check that nan in input array result in nan output + anan = np.array([[1,np.nan],[-1,1]]) + assert_equal(stats.ttest_1samp(anan, 0),([0, np.nan], [1,np.nan])) + finally: + np.seterr(**olderr) + + +def test_describe(): + x = np.vstack((np.ones((3,4)),2*np.ones((2,4)))) + nc, mmc = (5, ([ 1., 1., 1., 1.], [ 2., 2., 2., 2.])) + mc = np.array([ 1.4, 1.4, 1.4, 1.4]) + vc = np.array([ 0.3, 0.3, 0.3, 0.3]) + skc = [0.40824829046386357]*4 + kurtc = [-1.833333333333333]*4 + n, mm, m, v, sk, kurt = stats.describe(x) + assert_equal(n, nc) + assert_equal(mm, mmc) + assert_equal(m, mc) + assert_equal(v, vc) + assert_array_almost_equal(sk, skc, decimal=13) #not sure about precision + assert_array_almost_equal(kurt, kurtc, decimal=13) + n, mm, m, v, sk, kurt = stats.describe(x.T, axis=1) + assert_equal(n, nc) + assert_equal(mm, mmc) + assert_equal(m, mc) + assert_equal(v, vc) + assert_array_almost_equal(sk, skc, decimal=13) #not sure about precision + assert_array_almost_equal(kurt, kurtc, decimal=13) + +def test_normalitytests(): + # numbers verified with R: dagoTest in package fBasics + st_normal, st_skew, st_kurt = (3.92371918, 1.98078826, -0.01403734) + pv_normal, pv_skew, pv_kurt = (0.14059673, 0.04761502, 0.98880019) + x = np.array((-2,-1,0,1,2,3)*4)**2 + yield assert_array_almost_equal, stats.normaltest(x), (st_normal, pv_normal) + yield assert_array_almost_equal, stats.skewtest(x), (st_skew, pv_skew) + yield assert_array_almost_equal, stats.kurtosistest(x), (st_kurt, pv_kurt) + + +#class TestJarqueBera(TestCase): +# def test_jarque_bera_stats(self): +# np.random.seed(987654321) +# x = np.random.normal(0, 1, 100000) +# y = np.random.chisquare(10000, 100000) +# z = np.random.rayleigh(1, 100000) +# +# assert_(stats.jarque_bera(x)[1] > stats.jarque_bera(y)[1]) +# assert_(stats.jarque_bera(x)[1] > stats.jarque_bera(z)[1]) +# assert_(stats.jarque_bera(y)[1] > stats.jarque_bera(z)[1]) +# +# def test_jarque_bera_array_like(self): +# np.random.seed(987654321) +# x = np.random.normal(0, 1, 100000) +# +# JB1, p1 = stats.jarque_bera(list(x)) +# JB2, p2 = stats.jarque_bera(tuple(x)) +# JB3, p3 = stats.jarque_bera(x.reshape(2, 50000)) +# +# assert_(JB1 == JB2 == JB3) +# assert_(p1 == p2 == p3) +# +# def test_jarque_bera_size(self): +# assert_raises(ValueError, stats.jarque_bera, []) + + +def test_skewtest_too_few_samples(): + """Regression test for ticket #1492. + + skewtest requires at least 8 samples; 7 should raise a ValueError. + """ + x = np.arange(7.0) + assert_raises(ValueError, stats.skewtest, x) + +def test_kurtosistest_too_few_samples(): + """Regression test for ticket #1425. + + kurtosistest requires at least 5 samples; 4 should raise a ValueError. + """ + x = np.arange(4.0) + assert_raises(ValueError, stats.kurtosistest, x) + +def mannwhitneyu(): + x = np.array([ 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., + 1., 1., 1., 1., 1., 1., 1., 1., 2., 1., 1., 1., 1., 1., 1., 1., + 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., + 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., + 1., 1., 2., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., + 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., + 1., 1., 1., 1., 2., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., + 1., 1., 1., 1., 1., 2., 1., 1., 1., 1., 2., 1., 1., 2., 1., 1., + 2., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., + 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 2., 1., 1., 1., 1., + 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., + 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., + 2., 1., 1., 1., 1., 1., 1., 1., 1., 1., 2., 1., 1., 1., 1., 1., + 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 3., 1., 1., + 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., + 1., 1., 1., 1., 1., 1., 1.]) + + y = np.array([ 1., 1., 1., 1., 1., 1., 1., 2., 1., 2., 1., 1., 1., + 1., 2., 1., 1., 1., 2., 1., 1., 1., 1., 1., 2., 1., 1., 3., 1., + 1., 1., 1., 1., 1., 1., 1., 1., 1., 2., 1., 2., 1., 1., 1., 1., + 1., 1., 2., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., + 1., 1., 2., 1., 1., 1., 1., 1., 2., 2., 1., 1., 2., 1., 1., 2., + 1., 2., 1., 1., 1., 1., 2., 2., 1., 1., 1., 1., 1., 1., 1., 1., + 1., 1., 1., 1., 1., 1., 2., 1., 1., 1., 1., 1., 2., 2., 2., 1., + 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., + 1., 2., 1., 1., 2., 1., 1., 1., 1., 2., 1., 1., 1., 1., 1., 1., + 1., 1., 1., 1., 1., 1., 2., 1., 1., 1., 2., 1., 1., 1., 1., 1., + 1.]) + #p-value verified with matlab and R to 5 significant digits + assert_array_almost_equal(stats.stats.mannwhitneyu(x,y), + (16980.5, 2.8214327656317373e-005), decimal=12) + + + +def test_pointbiserial(): + # copied from mstats tests removing nans + x = [1,0,1,1,1,1,0,1,0,0,0,1,1,0,0,0,1,1,1,0,0,0,0,0,0,0,0,1,0, + 0,0,0,0,1] + y = [14.8,13.8,12.4,10.1,7.1,6.1,5.8,4.6,4.3,3.5,3.3,3.2,3.0, + 2.8,2.8,2.5,2.4,2.3,2.1,1.7,1.7,1.5,1.3,1.3,1.2,1.2,1.1, + 0.8,0.7,0.6,0.5,0.2,0.2,0.1] + assert_almost_equal(stats.pointbiserialr(x, y)[0], 0.36149, 5) + + +def test_obrientransform(): + #this is a regression test to check np.var replacement + #I didn't separately verigy the numbers + x1 = np.arange(5) + result = np.array( + [[ 5.41666667, 1.04166667, -0.41666667, 1.04166667, 5.41666667], + [ 21.66666667, 4.16666667, -1.66666667, 4.16666667, 21.66666667]]) + assert_array_almost_equal(stats.obrientransform(x1, 2*x1), result, decimal=8) + + +class HarMeanTestCase: + def test_1dlist(self): + ''' Test a 1d list''' + a=[10, 20, 30, 40, 50, 60, 70, 80, 90, 100] + b = 34.1417152147 + self.do(a, b) + def test_1darray(self): + ''' Test a 1d array''' + a=np.array([10, 20, 30, 40, 50, 60, 70, 80, 90, 100]) + b = 34.1417152147 + self.do(a, b) + def test_1dma(self): + ''' Test a 1d masked array''' + a=np.ma.array([10, 20, 30, 40, 50, 60, 70, 80, 90, 100]) + b = 34.1417152147 + self.do(a, b) + def test_1dmavalue(self): + ''' Test a 1d masked array with a masked value''' + a=np.ma.array([10, 20, 30, 40, 50, 60, 70, 80, 90, 100], + mask=[0,0,0,0,0,0,0,0,0,1]) + b = 31.8137186141 + self.do(a, b) + + # Note the next tests use axis=None as default, not axis=0 + def test_2dlist(self): + ''' Test a 2d list''' + a=[[10, 20, 30, 40], [50, 60, 70, 80], [90, 100, 110, 120]] + b = 38.6696271841 + self.do(a, b) + def test_2darray(self): + ''' Test a 2d array''' + a=[[10, 20, 30, 40], [50, 60, 70, 80], [90, 100, 110, 120]] + b = 38.6696271841 + self.do(np.array(a), b) + def test_2dma(self): + ''' Test a 2d masked array''' + a=[[10, 20, 30, 40], [50, 60, 70, 80], [90, 100, 110, 120]] + b = 38.6696271841 + self.do(np.ma.array(a), b) + def test_2daxis0(self): + ''' Test a 2d list with axis=0''' + a=[[10, 20, 30, 40], [50, 60, 70, 80], [90, 100, 110, 120]] + b = np.array([ 22.88135593, 39.13043478, 52.90076336, 65.45454545]) + self.do(a, b, axis=0) + def test_2daxis1(self): + ''' Test a 2d list with axis=1''' + a=[[10, 20, 30, 40], [50, 60, 70, 80], [90, 100, 110, 120]] + b = np.array([ 19.2 , 63.03939962, 103.80078637]) + self.do(a, b, axis=1) + def test_2dmatrixdaxis0(self): + ''' Test a 2d list with axis=0''' + a=[[10, 20, 30, 40], [50, 60, 70, 80], [90, 100, 110, 120]] + b = np.matrix([[ 22.88135593, 39.13043478, 52.90076336, 65.45454545]]) + self.do(np.matrix(a), b, axis=0) + def test_2dmatrixaxis1(self): + ''' Test a 2d list with axis=1''' + a=[[10, 20, 30, 40], [50, 60, 70, 80], [90, 100, 110, 120]] + b = np.matrix([[ 19.2 , 63.03939962, 103.80078637]]).T + self.do(np.matrix(a), b, axis=1) +## def test_dtype(self): +## ''' Test a 1d list with a new dtype''' +## a=[10, 20, 30, 40, 50, 60, 70, 80, 90, 100] +## b = 34.1417152147 +## self.do(a, b, dtype=np.float128) # does not work on Win32 + +class TestHarMean(HarMeanTestCase, TestCase): + def do(self, a, b, axis=None, dtype=None): + x = stats.hmean(a, axis=axis, dtype=dtype) + assert_almost_equal(b, x) + assert_equal(x.dtype, dtype) + +class GeoMeanTestCase: + def test_1dlist(self): + ''' Test a 1d list''' + a=[10, 20, 30, 40, 50, 60, 70, 80, 90, 100] + b = 45.2872868812 + self.do(a, b) + def test_1darray(self): + ''' Test a 1d array''' + a=np.array([10, 20, 30, 40, 50, 60, 70, 80, 90, 100]) + b = 45.2872868812 + self.do(a, b) + def test_1dma(self): + ''' Test a 1d masked array''' + a=np.ma.array([10, 20, 30, 40, 50, 60, 70, 80, 90, 100]) + b = 45.2872868812 + self.do(a, b) + def test_1dmavalue(self): + ''' Test a 1d masked array with a masked value''' + a=np.ma.array([10, 20, 30, 40, 50, 60, 70, 80, 90, 100], mask=[0,0,0,0,0,0,0,0,0,1]) + b = 41.4716627439 + self.do(a, b) + + # Note the next tests use axis=None as default, not axis=0 + def test_2dlist(self): + ''' Test a 2d list''' + a=[[10, 20, 30, 40], [50, 60, 70, 80], [90, 100, 110, 120]] + b = 52.8885199 + self.do(a, b) + def test_2darray(self): + ''' Test a 2d array''' + a=[[10, 20, 30, 40], [50, 60, 70, 80], [90, 100, 110, 120]] + b = 52.8885199 + self.do(np.array(a), b) + def test_2dma(self): + ''' Test a 2d masked array''' + a=[[10, 20, 30, 40], [50, 60, 70, 80], [90, 100, 110, 120]] + b = 52.8885199 + self.do(np.ma.array(a), b) + def test_2daxis0(self): + ''' Test a 2d list with axis=0''' + a=[[10, 20, 30, 40], [50, 60, 70, 80], [90, 100, 110, 120]] + b = np.array([35.56893304, 49.32424149, 61.3579244 , 72.68482371]) + self.do(a, b, axis=0) + def test_2daxis1(self): + ''' Test a 2d list with axis=1''' + a=[[10, 20, 30, 40], [50, 60, 70, 80], [90, 100, 110, 120]] + b = np.array([ 22.13363839, 64.02171746, 104.40086817]) + self.do(a, b, axis=1) + def test_2dmatrixdaxis0(self): + ''' Test a 2d list with axis=0''' + a=[[10, 20, 30, 40], [50, 60, 70, 80], [90, 100, 110, 120]] + b = np.matrix([[35.56893304, 49.32424149, 61.3579244 , 72.68482371]]) + self.do(np.matrix(a), b, axis=0) + def test_2dmatrixaxis1(self): + ''' Test a 2d list with axis=1''' + a=[[10, 20, 30, 40], [50, 60, 70, 80], [90, 100, 110, 120]] + b = np.matrix([[ 22.13363839, 64.02171746, 104.40086817]]).T + self.do(np.matrix(a), b, axis=1) +## def test_dtype(self): +## ''' Test a 1d list with a new dtype''' +## a=[10, 20, 30, 40, 50, 60, 70, 80, 90, 100] +## b = 45.2872868812 +## self.do(a, b, dtype=np.float128) # does not exist on win32 + def test_1dlist0(self): + ''' Test a 1d list with zero element''' + a=[10, 20, 30, 40, 50, 60, 70, 80, 90, 0] + b = 0.0 # due to exp(-inf)=0 + olderr = np.seterr(all='ignore') + try: + self.do(a, b) + finally: + np.seterr(**olderr) + + def test_1darray0(self): + ''' Test a 1d array with zero element''' + a=np.array([10, 20, 30, 40, 50, 60, 70, 80, 90, 0]) + b = 0.0 # due to exp(-inf)=0 + olderr = np.seterr(all='ignore') + try: + self.do(a, b) + finally: + np.seterr(**olderr) + + def test_1dma0(self): + ''' Test a 1d masked array with zero element''' + a=np.ma.array([10, 20, 30, 40, 50, 60, 70, 80, 90, 0]) + b = 41.4716627439 + olderr = np.seterr(all='ignore') + try: + self.do(a, b) + finally: + np.seterr(**olderr) + + def test_1dmainf(self): + ''' Test a 1d masked array with negative element''' + a=np.ma.array([10, 20, 30, 40, 50, 60, 70, 80, 90, -1]) + b = 41.4716627439 + olderr = np.seterr(all='ignore') + try: + self.do(a, b) + finally: + np.seterr(**olderr) + +class TestGeoMean(GeoMeanTestCase, TestCase): + def do(self, a, b, axis=None, dtype=None): + #Note this doesn't test when axis is not specified + x = stats.gmean(a, axis=axis, dtype=dtype) + assert_almost_equal(b, x) + assert_equal(x.dtype, dtype) + + +def test_binomtest(): + # precision tests compared to R for ticket:986 + pp = np.concatenate(( np.linspace(0.1,0.2,5), np.linspace(0.45,0.65,5), + np.linspace(0.85,0.95,5))) + n = 501 + x = 450 + results = [0.0, 0.0, 1.0159969301994141e-304, + 2.9752418572150531e-275, 7.7668382922535275e-250, + 2.3381250925167094e-099, 7.8284591587323951e-081, + 9.9155947819961383e-065, 2.8729390725176308e-050, + 1.7175066298388421e-037, 0.0021070691951093692, + 0.12044570587262322, 0.88154763174802508, 0.027120993063129286, + 2.6102587134694721e-006] + + for p, res in zip(pp,results): + assert_approx_equal(stats.binom_test(x, n, p), res, + significant=12, err_msg='fail forp=%f'%p) + + assert_approx_equal(stats.binom_test(50,100,0.1), 5.8320387857343647e-024, + significant=12, err_msg='fail forp=%f'%p) + +class Test_Trim(object): + # test trim functions + def test_trim1(self): + a = np.arange(11) + assert_equal(stats.trim1(a, 0.1), np.arange(10)) + assert_equal(stats.trim1(a, 0.2), np.arange(9)) + assert_equal(stats.trim1(a, 0.2, tail='left'), np.arange(2,11)) + assert_equal(stats.trim1(a, 3/11., tail='left'), np.arange(3,11)) + + def test_trimboth(self): + a = np.arange(11) + assert_equal(stats.trimboth(a, 3/11.), np.arange(3,8)) + assert_equal(stats.trimboth(a, 0.2), np.array([2, 3, 4, 5, 6, 7, 8])) + assert_equal(stats.trimboth(np.arange(24).reshape(6,4), 0.2), + np.arange(4,20).reshape(4,4)) + assert_equal(stats.trimboth(np.arange(24).reshape(4,6).T, 2/6.), + np.array([[ 2, 8, 14, 20],[ 3, 9, 15, 21]])) + assert_raises(ValueError, stats.trimboth, + np.arange(24).reshape(4,6).T, 4/6.) + + def test_trim_mean(self): + assert_equal(stats.trim_mean(np.arange(24).reshape(4,6).T, 2/6.), + np.array([ 2.5, 8.5, 14.5, 20.5])) + assert_equal(stats.trim_mean(np.arange(24).reshape(4,6), 2/6.), + np.array([ 9., 10., 11., 12., 13., 14.])) + assert_equal(stats.trim_mean(np.arange(24), 2/6.), 11.5) + assert_equal(stats.trim_mean([5,4,3,1,2,0], 2/6.), 2.5) + + +class TestSigamClip(object): + def test_sigmaclip1(self): + a = np.concatenate((np.linspace(9.5,10.5,31),np.linspace(0,20,5))) + fact = 4 #default + c, low, upp = stats.sigmaclip(a) + assert_(c.min()>low) + assert_(c.max()low) + assert_(c.max()low) + assert_(c.max()>> np = JITImport('numpy') - >>> np.exp(0)==1.0 - True - ''' -def test_bunch(): - ''' - >>> d = Bunch(test1=1,test2=3) - >>> d.test1; d.test2 - 1 - 3 - ''' -def test_dotdict(): - ''' - >>> d = DotDict(test1=1,test2=3) - >>> d.test1; d.test2 - 1 - 3 - ''' - -def test_detrendma(): - ''' - >>> x = linspace(0,1,200) - >>> y = exp(x)+0.1*cos(20*2*pi*x) - >>> y0 = detrendma(y,20); tr = y-y0 - >>> y0,tr - (array([ -1.05815186e-02, -2.48280355e-02, -7.01800760e-02, - -1.27193089e-01, -1.71915213e-01, -1.85125121e-01, - -1.59745361e-01, -1.03571981e-01, -3.62676515e-02, - 1.82219951e-02, 4.09039083e-02, 2.50630186e-02, - -2.11478040e-02, -7.78521440e-02, -1.21116040e-01, - -1.32178923e-01, -1.04689244e-01, -4.71541301e-02, - 2.03417510e-02, 7.38826137e-02, 8.95349902e-02, - 6.68738432e-02, 1.46828486e-02, -4.68648556e-02, - -9.39871606e-02, -1.08465407e-01, -8.46710629e-02, - -3.17365657e-02, 2.99669288e-02, 7.66864134e-02, - 9.04482283e-02, 6.59902473e-02, 1.27914062e-02, - -4.85841870e-02, -9.44185349e-02, -1.06987444e-01, - -8.13964951e-02, -2.74687460e-02, 3.40438793e-02, - 7.94643163e-02, 9.13222681e-02, 6.50922520e-02, - 1.09390148e-02, -5.02028639e-02, -9.47031411e-02, - -1.05349757e-01, -7.79872833e-02, -2.31196073e-02, - 3.81412653e-02, 8.22178144e-02, 9.21605209e-02, - 6.41850565e-02, 9.13184690e-03, -5.17149253e-02, - -9.48363260e-02, -1.03549587e-01, -7.44424124e-02, - -1.86890490e-02, 4.22594607e-02, 8.49486437e-02, - 9.29666543e-02, 6.32740911e-02, 7.37625254e-03, - -5.31142920e-02, -9.48133620e-02, -1.01584110e-01, - -7.07607748e-02, -1.41768231e-02, 4.63990484e-02, - 8.76587937e-02, 9.37446001e-02, 6.23650231e-02, - 5.67876495e-03, -5.43947621e-02, -9.46294406e-02, - -9.94504301e-02, -6.69411601e-02, -9.58252265e-03, - 5.05608316e-02, 9.03505172e-02, 9.44985623e-02, - 6.14637631e-02, 4.04610591e-03, -5.55500040e-02, - -9.42796647e-02, -9.71455674e-02, -6.29822440e-02, - -4.90556961e-03, 5.47458452e-02, 9.30263409e-02, - 9.52330253e-02, 6.05764719e-02, 2.48519180e-03, - -5.65735506e-02, -9.37590405e-02, -9.46664506e-02, - -5.88825766e-02, -1.45202622e-04, 5.89553685e-02, - 9.56890756e-02, 9.59527629e-02, 5.97095676e-02, - 1.00314001e-03, -5.74587921e-02, -9.30624694e-02, - -9.20099048e-02, -5.46405701e-02, 4.69953603e-03, - 6.31909369e-02, 9.83418277e-02, 9.66628470e-02, - 5.88697331e-02, -3.92724035e-04, -5.81989687e-02, - -9.21847386e-02, -8.91726414e-02, -5.02544862e-02, - 9.62981387e-03, 6.74543554e-02, 1.00988010e-01, - 9.73686580e-02, 5.80639242e-02, -1.69485946e-03, - -5.87871620e-02, -9.11205115e-02, -8.61512458e-02, - -4.57224228e-02, 1.46470222e-02, 7.17477118e-02, - 1.03631355e-01, 9.80758938e-02, 5.72993780e-02, - -2.89550192e-03, -5.92162868e-02, -8.98643173e-02, - -8.29421650e-02, -4.10422999e-02, 1.97527907e-02, - 7.60733908e-02, 1.06275925e-01, 9.87905812e-02, - 5.65836223e-02, -3.98665495e-03, -5.94790815e-02, - -8.84105398e-02, -7.95416952e-02, -3.62118451e-02, - 2.49490024e-02, 8.04340881e-02, 1.08926128e-01, - 9.95190863e-02, 5.59244846e-02, -4.96008086e-03, - -5.95680980e-02, -8.67534061e-02, -7.59459673e-02, - -3.12285782e-02, 3.02378095e-02, 8.48328258e-02, - 1.11586726e-01, 1.00268126e-01, 5.53301029e-02, - -5.80729079e-03, -5.94756912e-02, -8.48869734e-02, - -7.21509330e-02, -2.60897955e-02, 3.56216499e-02, - 8.92729678e-02, 1.14262857e-01, 1.01044781e-01, - 5.48089359e-02, -6.51953427e-03, -5.91940075e-02, - -8.28051165e-02, -6.81523491e-02, -2.07925530e-02, - 4.11032641e-02, 9.37582360e-02, 1.16960041e-01, - 1.13824241e-01, 7.82451609e-02, 2.87461256e-02, - -1.07566250e-02, -2.01779675e-02, 8.98967999e-03, - 7.03952281e-02, 1.45278564e-01, 2.09706186e-01, - 2.43802139e-01, 2.39414013e-01, 2.03257341e-01, - 1.54325635e-01, 1.16564992e-01, 1.09638547e-01, - 1.41342814e-01, 2.04600808e-01, 2.80191671e-01, - 3.44164010e-01, 3.77073744e-01]), array([ 1.11058152, 1.11058152, 1.11058152, 1.11058152, 1.11058152, - 1.11058152, 1.11058152, 1.11058152, 1.11058152, 1.11058152, - 1.11058152, 1.11058152, 1.11058152, 1.11058152, 1.11058152, - 1.11058152, 1.11058152, 1.11058152, 1.11058152, 1.11058152, - 1.11599212, 1.12125245, 1.12643866, 1.13166607, 1.13704477, - 1.14263723, 1.14843422, 1.15435845, 1.16029443, 1.16613308, - 1.17181383, 1.17734804, 1.18281471, 1.18833001, 1.19400259, - 1.19989168, 1.20598434, 1.21220048, 1.21842384, 1.22454684, - 1.23051218, 1.23633498, 1.24209697, 1.24791509, 1.25389641, - 1.26009689, 1.26649987, 1.27302256, 1.27954802, 1.28597031, - 1.29223546, 1.29836228, 1.30443522, 1.31057183, 1.31687751, - 1.32340488, 1.3301336 , 1.33697825, 1.34382132, 1.35055864, - 1.35713958, 1.36358668, 1.36998697, 1.37645853, 1.38310497, - 1.38997553, 1.39704621, 1.40422902, 1.41140604, 1.41847493, - 1.4253885 , 1.43217295, 1.43891784, 1.44574164, 1.45274607, - 1.45997696, 1.46740665, 1.47494469, 1.48247285, 1.48989073, - 1.49715462, 1.50429437, 1.51140198, 1.51859618, 1.52597672, - 1.53358594, 1.54139257, 1.5493038 , 1.55720119, 1.56498641, - 1.57261924, 1.58013316, 1.58762252, 1.5952062 , 1.60298187, - 1.61098836, 1.6191908 , 1.62749412, 1.63577979, 1.64395163, - 1.65197298, 1.65988092, 1.66777202, 1.67576523, 1.68395602, - 1.69237968, 1.70099778, 1.70971307, 1.71840707, 1.72698583, - 1.73541631, 1.74373911, 1.75205298, 1.76047677, 1.76910369, - 1.77796544, 1.78702008, 1.79616827, 1.80529169, 1.81429875, - 1.82316 , 1.83191959, 1.84067831, 1.84955481, 1.85863994, - 1.86796178, 1.87747491, 1.88707803, 1.89665308, 1.9061109 , - 1.91542572, 1.92464514, 1.9338719 , 1.94322436, 1.9527909 , - 1.96259596, 1.97259069, 1.9826719 , 1.99272195, 2.00265419, - 2.01244653, 2.02215 , 2.0318692 , 2.04172204, 2.05179437, - 2.06210696, 2.07260759, 2.08319129, 2.09374092, 2.10417247, - 2.11446752, 2.12468051, 2.13491776, 2.14529665, 2.1559004 , - 2.16674609, 2.17777817, 2.18889002, 2.19996511, 2.21092214, - 2.22174641, 2.23249567, 2.24327791, 2.25420982, 2.26537192, - 2.2767776 , 2.28836802, 2.30003501, 2.3116628 , 2.32317284, - 2.33455419, 2.34586786, 2.35722337, 2.36873665, 2.38048542, - 2.39247934, 2.4046564 , 2.41690694, 2.42911606, 2.44120808, - 2.44120808, 2.44120808, 2.44120808, 2.44120808, 2.44120808, - 2.44120808, 2.44120808, 2.44120808, 2.44120808, 2.44120808, - 2.44120808, 2.44120808, 2.44120808, 2.44120808, 2.44120808, - 2.44120808, 2.44120808, 2.44120808, 2.44120808, 2.44120808])) - ''' - -def test_findcross_and_ecross(): - ''' - >>> findcross([0, 0, 1, -1, 1],0) - array([1, 2, 3]) - >>> findcross([0, 1, -1, 1],0) - array([0, 1, 2]) - - >>> t = linspace(0,7*pi,250) - >>> x = sin(t) - >>> ind = findcross(x,0.75) - >>> ind - array([ 9, 25, 80, 97, 151, 168, 223, 239]) - >>> t0 = ecross(t,x,ind,0.75) - >>> t0 - array([ 0.84910514, 2.2933879 , 7.13205663, 8.57630119, - 13.41484739, 14.85909194, 19.69776067, 21.14204343]) - - ''' -def test_findextrema(): - ''' - >>> t = linspace(0,7*pi,250) - >>> x = sin(t) - >>> ind = findextrema(x) - >>> ind - array([ 18, 53, 89, 125, 160, 196, 231]) - ''' - -def test_findrfc(): - ''' - - >>> t = linspace(0,7*pi,250) - >>> x = sin(t)+0.1*sin(50*t) - >>> ind = findextrema(x) - >>> ind - array([ 1, 3, 4, 6, 7, 9, 11, 13, 14, 16, 18, 19, 21, - 23, 25, 26, 28, 29, 31, 33, 35, 36, 38, 39, 41, 43, - 45, 46, 48, 50, 51, 53, 55, 56, 58, 60, 61, 63, 65, - 67, 68, 70, 71, 73, 75, 77, 78, 80, 81, 83, 85, 87, - 88, 90, 92, 93, 95, 97, 99, 100, 102, 103, 105, 107, 109, - 110, 112, 113, 115, 117, 119, 120, 122, 124, 125, 127, 129, 131, - 132, 134, 135, 137, 139, 141, 142, 144, 145, 147, 149, 151, 152, - 154, 156, 157, 159, 161, 162, 164, 166, 167, 169, 171, 173, 174, - 176, 177, 179, 181, 183, 184, 186, 187, 189, 191, 193, 194, 196, - 198, 199, 201, 203, 205, 206, 208, 209, 211, 213, 215, 216, 218, - 219, 221, 223, 225, 226, 228, 230, 231, 233, 235, 237, 238, 240, - 241, 243, 245, 247, 248]) - >>> ti, tp = t[ind], x[ind] - >>> ind1 = findrfc(tp,0.3) - >>> ind1 - array([ 0, 9, 32, 53, 74, 95, 116, 137]) - >>> tp[ind1] - array([-0.00743352, 1.08753972, -1.07206545, 1.09550837, -1.07940458, - 1.07849396, -1.0995006 , 1.08094452]) - ''' - -def test_rfcfilter(): - ''' - # 1. Filtered signal y is the turning points of x. - >>> x = sea() - >>> y = rfcfilter(x[:,1], h=0, method=1) - >>> y[0:5] - array([-1.2004945 , 0.83950546, -0.09049454, -0.02049454, -0.09049454]) - - # 2. This removes all rainflow cycles with range less than 0.5. - >>> y1 = rfcfilter(x[:,1], h=0.5) - >>> y1[0:5] - array([-1.2004945 , 0.83950546, -0.43049454, 0.34950546, -0.51049454]) - - >>> t = linspace(0,7*pi,250) - >>> x = sin(t)+0.1*sin(50*t) - >>> ind = findextrema(x) - >>> ind - array([ 1, 3, 4, 6, 7, 9, 11, 13, 14, 16, 18, 19, 21, - 23, 25, 26, 28, 29, 31, 33, 35, 36, 38, 39, 41, 43, - 45, 46, 48, 50, 51, 53, 55, 56, 58, 60, 61, 63, 65, - 67, 68, 70, 71, 73, 75, 77, 78, 80, 81, 83, 85, 87, - 88, 90, 92, 93, 95, 97, 99, 100, 102, 103, 105, 107, 109, - 110, 112, 113, 115, 117, 119, 120, 122, 124, 125, 127, 129, 131, - 132, 134, 135, 137, 139, 141, 142, 144, 145, 147, 149, 151, 152, - 154, 156, 157, 159, 161, 162, 164, 166, 167, 169, 171, 173, 174, - 176, 177, 179, 181, 183, 184, 186, 187, 189, 191, 193, 194, 196, - 198, 199, 201, 203, 205, 206, 208, 209, 211, 213, 215, 216, 218, - 219, 221, 223, 225, 226, 228, 230, 231, 233, 235, 237, 238, 240, - 241, 243, 245, 247, 248]) - >>> ti, tp = t[ind], x[ind] - >>> tp03 = rfcfilter(tp,0.3) - >>> tp03 - array([-0.00743352, 1.08753972, -1.07206545, 1.09550837, -1.07940458, - 1.07849396, -1.0995006 , 1.08094452, 0.11983423]) - ''' -def test_findtp(): - ''' - >>> import numpy as np - >>> x = sea() - >>> x1 = x[0:200,:] - >>> itp = findtp(x1[:,1],0,'Mw') - >>> itph = findtp(x1[:,1],0.3,'Mw') - >>> itp - array([ 11, 21, 22, 24, 26, 28, 31, 39, 43, 45, 47, 51, 56, - 64, 70, 78, 82, 84, 89, 94, 101, 108, 119, 131, 141, 148, - 149, 150, 159, 173, 184, 190, 199]) - >>> itph - array([ 11, 28, 31, 39, 47, 51, 56, 64, 70, 78, 89, 94, 101, - 108, 119, 131, 141, 148, 159, 173, 184, 190, 199]) - ''' -def test_findtc(): - ''' - >>> x = sea() - >>> x1 = x[0:200,:] - >>> itc, iv = findtc(x1[:,1],0,'dw') - >>> itc - array([ 28, 31, 39, 56, 64, 69, 78, 82, 83, 89, 94, 101, 108, - 119, 131, 140, 148, 159, 173, 184]) - >>> iv - array([ 19, 29, 34, 53, 60, 67, 76, 81, 82, 84, 90, 99, 103, - 112, 127, 137, 143, 154, 166, 180, 185]) - ''' - -def test_findoutliers(): - ''' - >>> xx = sea() - >>> dt = diff(xx[:2,0]) - >>> dcrit = 5*dt - >>> ddcrit = 9.81/2*dt*dt - >>> zcrit = 0 - >>> [inds, indg] = findoutliers(xx[:,1],zcrit,dcrit,ddcrit,verbose=True) - Found 0 spurious positive jumps of Dx - Found 0 spurious negative jumps of Dx - Found 37 spurious positive jumps of D^2x - Found 200 spurious negative jumps of D^2x - Found 244 consecutive equal values - Found the total of 1152 spurious points - >>> inds - array([ 6, 7, 8, ..., 9509, 9510, 9511]) - >>> indg - array([ 0, 1, 2, ..., 9521, 9522, 9523]) - ''' -def test_common_shape(): - ''' - >>> import numpy as np - >>> A = np.ones((4,1)) - >>> B = 2 - >>> C = np.ones((1,5))*5 - >>> common_shape(A,B,C) - (4, 5) - >>> common_shape(A,B,C,shape=(3,4,1)) - (3, 4, 5) - >>> A = np.ones((4,1)) - >>> B = 2 - >>> C = np.ones((1,5))*5 - >>> common_shape(A,B,C) - (4, 5) - >>> common_shape(A,B,C,shape=(3,4,1)) - (3, 4, 5) - ''' -def test_argsreduce(): - ''' - >>> import numpy as np - >>> rand = np.random.random_sample - >>> A = linspace(0,19,20).reshape((4,5)) - >>> B = 2 - >>> C = range(5) - >>> cond = np.ones(A.shape) - >>> [A1,B1,C1] = argsreduce(cond,A,B,C) - >>> B1.shape - (20,) - >>> cond[2,:] = 0 - >>> [A2,B2,C2] = argsreduce(cond,A,B,C) - >>> B2.shape - (15,) - >>> A2;B2;C2 - array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 15., - 16., 17., 18., 19.]) - array([2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]) - array([0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4]) - ''' -def test_stirlerr(): - ''' - >>> stirlerr(range(5)) - array([ inf, 0.08106147, 0.0413407 , 0.02767793, 0.02079067]) - ''' -def test_getshipchar(): - ''' - >>> sc = getshipchar(10,'service_speed') - >>> names = ['beam', 'beamSTD', 'draught', - ... 'draughtSTD', 'length', 'lengthSTD', - ... 'max_deadweight', 'max_deadweightSTD', 'propeller_diameter', - ... 'propeller_diameterSTD', 'service_speed', 'service_speedSTD'] - >>> for name in names: print( '%s : %g' % (name, sc[name])) - beam : 29 - beamSTD : 2.9 - draught : 9.6 - draughtSTD : 2.112 - length : 216 - lengthSTD : 2.01131 - max_deadweight : 30969 - max_deadweightSTD : 3096.9 - propeller_diameter : 6.76117 - propeller_diameterSTD : 0.20267 - service_speed : 10 - service_speedSTD : 0 - ''' -def test_betaloge(): - ''' - >>> betaloge(3, arange(4)) - array([ inf, -1.09861229, -2.48490665, -3.40119738]) - ''' -def test_gravity(): - ''' - >>> phi = linspace(0,45,5) - >>> gravity(phi) - array([ 9.78049 , 9.78245014, 9.78803583, 9.79640552, 9.80629387]) - ''' -def test_nextpow2(): - ''' - >>> nextpow2(10) - 4 - >>> nextpow2(np.arange(5)) - 3 - ''' - -def test_discretize(): - ''' - >>> x, y = discretize(np.cos,0,np.pi) - >>> x; y - array([ 0. , 0.19634954, 0.39269908, 0.58904862, 0.78539816, - 0.9817477 , 1.17809725, 1.37444679, 1.57079633, 1.76714587, - 1.96349541, 2.15984495, 2.35619449, 2.55254403, 2.74889357, - 2.94524311, 3.14159265]) - array([ 1.00000000e+00, 9.80785280e-01, 9.23879533e-01, - 8.31469612e-01, 7.07106781e-01, 5.55570233e-01, - 3.82683432e-01, 1.95090322e-01, 6.12323400e-17, - -1.95090322e-01, -3.82683432e-01, -5.55570233e-01, - -7.07106781e-01, -8.31469612e-01, -9.23879533e-01, - -9.80785280e-01, -1.00000000e+00]) - ''' -def test_discretize_adaptive(): - ''' - >>> x, y = discretize(np.cos,0,np.pi, method='adaptive') - >>> x; y - array([ 0. , 0.19634954, 0.39269908, 0.58904862, 0.78539816, - 0.9817477 , 1.17809725, 1.37444679, 1.57079633, 1.76714587, - 1.96349541, 2.15984495, 2.35619449, 2.55254403, 2.74889357, - 2.94524311, 3.14159265]) - array([ 1.00000000e+00, 9.80785280e-01, 9.23879533e-01, - 8.31469612e-01, 7.07106781e-01, 5.55570233e-01, - 3.82683432e-01, 1.95090322e-01, 6.12323400e-17, - -1.95090322e-01, -3.82683432e-01, -5.55570233e-01, - -7.07106781e-01, -8.31469612e-01, -9.23879533e-01, - -9.80785280e-01, -1.00000000e+00]) - ''' -def test_pol2cart_n_cart2pol(): - ''' - >>> r = 5 - >>> t = linspace(0,pi,20) - >>> x, y = polar2cart(t,r) - >>> x; y - array([ 5. , 4.93180652, 4.72908621, 4.39736876, 3.94570255, - 3.38640786, 2.73474079, 2.00847712, 1.22742744, 0.41289673, - -0.41289673, -1.22742744, -2.00847712, -2.73474079, -3.38640786, - -3.94570255, -4.39736876, -4.72908621, -4.93180652, -5. ]) - array([ 0.00000000e+00, 8.22972951e-01, 1.62349735e+00, - 2.37973697e+00, 3.07106356e+00, 3.67861955e+00, - 4.18583239e+00, 4.57886663e+00, 4.84700133e+00, - 4.98292247e+00, 4.98292247e+00, 4.84700133e+00, - 4.57886663e+00, 4.18583239e+00, 3.67861955e+00, - 3.07106356e+00, 2.37973697e+00, 1.62349735e+00, - 8.22972951e-01, 6.12323400e-16]) - >>> ti, ri = cart2polar(x,y) - >>> ti;ri - array([ 0. , 0.16534698, 0.33069396, 0.49604095, 0.66138793, - 0.82673491, 0.99208189, 1.15742887, 1.32277585, 1.48812284, - 1.65346982, 1.8188168 , 1.98416378, 2.14951076, 2.31485774, - 2.48020473, 2.64555171, 2.81089869, 2.97624567, 3.14159265]) - array([ 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., - 5., 5., 5., 5., 5., 5., 5.]) - ''' -def test_meshgrid(): - ''' - >>> x = np.linspace(0,1,3) # coordinates along x axis - >>> y = np.linspace(0,1,2) # coordinates along y axis - >>> xv, yv = meshgrid(x,y) # extend x and y for a 2D xy grid - >>> xv - array([[ 0. , 0.5, 1. ], - [ 0. , 0.5, 1. ]]) - >>> yv - array([[ 0., 0., 0.], - [ 1., 1., 1.]]) - >>> xv, yv = meshgrid(x,y, sparse=True) # make sparse output arrays - >>> xv - array([[ 0. , 0.5, 1. ]]) - >>> yv - array([[ 0.], - [ 1.]]) - - >>> meshgrid(x,y,sparse=True,indexing='ij') # change to matrix indexing - [array([[ 0. ], - [ 0.5], - [ 1. ]]), array([[ 0., 1.]])] - >>> meshgrid(x,y,indexing='ij') - [array([[ 0. , 0. ], - [ 0.5, 0.5], - [ 1. , 1. ]]), array([[ 0., 1.], - [ 0., 1.], - [ 0., 1.]])] - - >>> meshgrid(0,1,5) # just a 3D point - [array([[[0]]]), array([[[1]]]), array([[[5]]])] - >>> map(np.squeeze,meshgrid(0,1,5)) # just a 3D point - [array(0), array(1), array(5)] - >>> meshgrid(3) - array([3]) - >>> meshgrid(y) # 1D grid y is just returned - array([ 0., 1.]) - - `meshgrid` is very useful to evaluate functions on a grid. - - >>> x = np.arange(-5, 5, 0.1) - >>> y = np.arange(-5, 5, 0.1) - >>> xx, yy = meshgrid(x, y, sparse=True) - >>> z = np.sin(xx**2+yy**2)/(xx**2+yy**2) - ''' -def test_tranproc(): - ''' - >>> import wafo.transform.models as wtm - >>> tr = wtm.TrHermite() - >>> x = linspace(-5,5,501) - >>> g = tr(x) - >>> y0, y1 = tranproc(x, g, range(5), ones(5)) - >>> y0;y1 - array([ 0.02659612, 1.00115284, 1.92872532, 2.81453257, 3.66292878]) - array([ 1.00005295, 0.9501118 , 0.90589954, 0.86643821, 0.83096482]) - - ''' -if __name__=='__main__': - import doctest - doctest.testmod() +import numpy as np # @UnusedImport +#@UnusedImport +from numpy import cos, exp, linspace, pi, sin, diff, arange, ones +from numpy.random import randn # @UnusedImport +from wafo.data import sea # @UnusedImport +from wafo.misc import (JITImport, Bunch, detrendma, DotDict, findcross, ecross, findextrema, # @UnusedImport + #@UnusedImport + findrfc, rfcfilter, findtp, findtc, findoutliers, + common_shape, argsreduce, stirlerr, getshipchar, betaloge, + #@UnusedImport + #@UnusedImport + gravity, nextpow2, discretize, polar2cart, + cart2polar, meshgrid, tranproc) # @UnusedImport + + +def test_JITImport(): + ''' + >>> np = JITImport('numpy') + >>> np.exp(0)==1.0 + True + ''' + + +def test_bunch(): + ''' + >>> d = Bunch(test1=1,test2=3) + >>> d.test1; d.test2 + 1 + 3 + ''' + + +def test_dotdict(): + ''' + >>> d = DotDict(test1=1,test2=3) + >>> d.test1; d.test2 + 1 + 3 + ''' + + +def test_detrendma(): + ''' + >>> x = linspace(0,1,200) + >>> y = exp(x)+0.1*cos(20*2*pi*x) + >>> y0 = detrendma(y,20); tr = y-y0 + >>> y0,tr + (array([ -1.05815186e-02, -2.48280355e-02, -7.01800760e-02, + -1.27193089e-01, -1.71915213e-01, -1.85125121e-01, + -1.59745361e-01, -1.03571981e-01, -3.62676515e-02, + 1.82219951e-02, 4.09039083e-02, 2.50630186e-02, + -2.11478040e-02, -7.78521440e-02, -1.21116040e-01, + -1.32178923e-01, -1.04689244e-01, -4.71541301e-02, + 2.03417510e-02, 7.38826137e-02, 8.95349902e-02, + 6.68738432e-02, 1.46828486e-02, -4.68648556e-02, + -9.39871606e-02, -1.08465407e-01, -8.46710629e-02, + -3.17365657e-02, 2.99669288e-02, 7.66864134e-02, + 9.04482283e-02, 6.59902473e-02, 1.27914062e-02, + -4.85841870e-02, -9.44185349e-02, -1.06987444e-01, + -8.13964951e-02, -2.74687460e-02, 3.40438793e-02, + 7.94643163e-02, 9.13222681e-02, 6.50922520e-02, + 1.09390148e-02, -5.02028639e-02, -9.47031411e-02, + -1.05349757e-01, -7.79872833e-02, -2.31196073e-02, + 3.81412653e-02, 8.22178144e-02, 9.21605209e-02, + 6.41850565e-02, 9.13184690e-03, -5.17149253e-02, + -9.48363260e-02, -1.03549587e-01, -7.44424124e-02, + -1.86890490e-02, 4.22594607e-02, 8.49486437e-02, + 9.29666543e-02, 6.32740911e-02, 7.37625254e-03, + -5.31142920e-02, -9.48133620e-02, -1.01584110e-01, + -7.07607748e-02, -1.41768231e-02, 4.63990484e-02, + 8.76587937e-02, 9.37446001e-02, 6.23650231e-02, + 5.67876495e-03, -5.43947621e-02, -9.46294406e-02, + -9.94504301e-02, -6.69411601e-02, -9.58252265e-03, + 5.05608316e-02, 9.03505172e-02, 9.44985623e-02, + 6.14637631e-02, 4.04610591e-03, -5.55500040e-02, + -9.42796647e-02, -9.71455674e-02, -6.29822440e-02, + -4.90556961e-03, 5.47458452e-02, 9.30263409e-02, + 9.52330253e-02, 6.05764719e-02, 2.48519180e-03, + -5.65735506e-02, -9.37590405e-02, -9.46664506e-02, + -5.88825766e-02, -1.45202622e-04, 5.89553685e-02, + 9.56890756e-02, 9.59527629e-02, 5.97095676e-02, + 1.00314001e-03, -5.74587921e-02, -9.30624694e-02, + -9.20099048e-02, -5.46405701e-02, 4.69953603e-03, + 6.31909369e-02, 9.83418277e-02, 9.66628470e-02, + 5.88697331e-02, -3.92724035e-04, -5.81989687e-02, + -9.21847386e-02, -8.91726414e-02, -5.02544862e-02, + 9.62981387e-03, 6.74543554e-02, 1.00988010e-01, + 9.73686580e-02, 5.80639242e-02, -1.69485946e-03, + -5.87871620e-02, -9.11205115e-02, -8.61512458e-02, + -4.57224228e-02, 1.46470222e-02, 7.17477118e-02, + 1.03631355e-01, 9.80758938e-02, 5.72993780e-02, + -2.89550192e-03, -5.92162868e-02, -8.98643173e-02, + -8.29421650e-02, -4.10422999e-02, 1.97527907e-02, + 7.60733908e-02, 1.06275925e-01, 9.87905812e-02, + 5.65836223e-02, -3.98665495e-03, -5.94790815e-02, + -8.84105398e-02, -7.95416952e-02, -3.62118451e-02, + 2.49490024e-02, 8.04340881e-02, 1.08926128e-01, + 9.95190863e-02, 5.59244846e-02, -4.96008086e-03, + -5.95680980e-02, -8.67534061e-02, -7.59459673e-02, + -3.12285782e-02, 3.02378095e-02, 8.48328258e-02, + 1.11586726e-01, 1.00268126e-01, 5.53301029e-02, + -5.80729079e-03, -5.94756912e-02, -8.48869734e-02, + -7.21509330e-02, -2.60897955e-02, 3.56216499e-02, + 8.92729678e-02, 1.14262857e-01, 1.01044781e-01, + 5.48089359e-02, -6.51953427e-03, -5.91940075e-02, + -8.28051165e-02, -6.81523491e-02, -2.07925530e-02, + 4.11032641e-02, 9.37582360e-02, 1.16960041e-01, + 1.13824241e-01, 7.82451609e-02, 2.87461256e-02, + -1.07566250e-02, -2.01779675e-02, 8.98967999e-03, + 7.03952281e-02, 1.45278564e-01, 2.09706186e-01, + 2.43802139e-01, 2.39414013e-01, 2.03257341e-01, + 1.54325635e-01, 1.16564992e-01, 1.09638547e-01, + 1.41342814e-01, 2.04600808e-01, 2.80191671e-01, + 3.44164010e-01, 3.77073744e-01]), array([ + 1.11058152, 1.11058152, 1.11058152, 1.11058152, 1.11058152, + 1.11058152, 1.11058152, 1.11058152, 1.11058152, 1.11058152, + 1.11058152, 1.11058152, 1.11058152, 1.11058152, 1.11058152, + 1.11058152, 1.11058152, 1.11058152, 1.11058152, 1.11058152, + 1.11599212, 1.12125245, 1.12643866, 1.13166607, 1.13704477, + 1.14263723, 1.14843422, 1.15435845, 1.16029443, 1.16613308, + 1.17181383, 1.17734804, 1.18281471, 1.18833001, 1.19400259, + 1.19989168, 1.20598434, 1.21220048, 1.21842384, 1.22454684, + 1.23051218, 1.23633498, 1.24209697, 1.24791509, 1.25389641, + 1.26009689, 1.26649987, 1.27302256, 1.27954802, 1.28597031, + 1.29223546, 1.29836228, 1.30443522, 1.31057183, 1.31687751, + 1.32340488, 1.3301336 , 1.33697825, 1.34382132, 1.35055864, + 1.35713958, 1.36358668, 1.36998697, 1.37645853, 1.38310497, + 1.38997553, 1.39704621, 1.40422902, 1.41140604, 1.41847493, + 1.4253885 , 1.43217295, 1.43891784, 1.44574164, 1.45274607, + 1.45997696, 1.46740665, 1.47494469, 1.48247285, 1.48989073, + 1.49715462, 1.50429437, 1.51140198, 1.51859618, 1.52597672, + 1.53358594, 1.54139257, 1.5493038 , 1.55720119, 1.56498641, + 1.57261924, 1.58013316, 1.58762252, 1.5952062 , 1.60298187, + 1.61098836, 1.6191908 , 1.62749412, 1.63577979, 1.64395163, + 1.65197298, 1.65988092, 1.66777202, 1.67576523, 1.68395602, + 1.69237968, 1.70099778, 1.70971307, 1.71840707, 1.72698583, + 1.73541631, 1.74373911, 1.75205298, 1.76047677, 1.76910369, + 1.77796544, 1.78702008, 1.79616827, 1.80529169, 1.81429875, + 1.82316 , 1.83191959, 1.84067831, 1.84955481, 1.85863994, + 1.86796178, 1.87747491, 1.88707803, 1.89665308, 1.9061109 , + 1.91542572, 1.92464514, 1.9338719 , 1.94322436, 1.9527909 , + 1.96259596, 1.97259069, 1.9826719 , 1.99272195, 2.00265419, + 2.01244653, 2.02215 , 2.0318692 , 2.04172204, 2.05179437, + 2.06210696, 2.07260759, 2.08319129, 2.09374092, 2.10417247, + 2.11446752, 2.12468051, 2.13491776, 2.14529665, 2.1559004 , + 2.16674609, 2.17777817, 2.18889002, 2.19996511, 2.21092214, + 2.22174641, 2.23249567, 2.24327791, 2.25420982, 2.26537192, + 2.2767776 , 2.28836802, 2.30003501, 2.3116628 , 2.32317284, + 2.33455419, 2.34586786, 2.35722337, 2.36873665, 2.38048542, + 2.39247934, 2.4046564 , 2.41690694, 2.42911606, 2.44120808, + 2.44120808, 2.44120808, 2.44120808, 2.44120808, 2.44120808, + 2.44120808, 2.44120808, 2.44120808, 2.44120808, 2.44120808, + 2.44120808, 2.44120808, 2.44120808, 2.44120808, 2.44120808, + 2.44120808, 2.44120808, 2.44120808, 2.44120808, 2.44120808])) + ''' + + +def test_findcross_and_ecross(): + ''' + >>> findcross([0, 0, 1, -1, 1],0) + array([1, 2, 3]) + >>> findcross([0, 1, -1, 1],0) + array([0, 1, 2]) + + >>> t = linspace(0,7*pi,250) + >>> x = sin(t) + >>> ind = findcross(x,0.75) + >>> ind + array([ 9, 25, 80, 97, 151, 168, 223, 239]) + >>> t0 = ecross(t,x,ind,0.75) + >>> t0 + array([ 0.84910514, 2.2933879 , 7.13205663, 8.57630119, + 13.41484739, 14.85909194, 19.69776067, 21.14204343]) + ''' + + +def test_findextrema(): + ''' + >>> t = linspace(0,7*pi,250) + >>> x = sin(t) + >>> ind = findextrema(x) + >>> ind + array([ 18, 53, 89, 125, 160, 196, 231]) + ''' + + +def test_findrfc(): + ''' + >>> t = linspace(0,7*pi,250) + >>> x = sin(t)+0.1*sin(50*t) + >>> ind = findextrema(x) + >>> ind + array([ 1, 3, 4, 6, 7, 9, 11, 13, 14, 16, 18, 19, 21, + 23, 25, 26, 28, 29, 31, 33, 35, 36, 38, 39, 41, 43, + 45, 46, 48, 50, 51, 53, 55, 56, 58, 60, 61, 63, 65, + 67, 68, 70, 71, 73, 75, 77, 78, 80, 81, 83, 85, 87, + 88, 90, 92, 93, 95, 97, 99, 100, 102, 103, 105, 107, 109, + 110, 112, 113, 115, 117, 119, 120, 122, 124, 125, 127, 129, 131, + 132, 134, 135, 137, 139, 141, 142, 144, 145, 147, 149, 151, 152, + 154, 156, 157, 159, 161, 162, 164, 166, 167, 169, 171, 173, 174, + 176, 177, 179, 181, 183, 184, 186, 187, 189, 191, 193, 194, 196, + 198, 199, 201, 203, 205, 206, 208, 209, 211, 213, 215, 216, 218, + 219, 221, 223, 225, 226, 228, 230, 231, 233, 235, 237, 238, 240, + 241, 243, 245, 247, 248]) + >>> ti, tp = t[ind], x[ind] + >>> ind1 = findrfc(tp,0.3) + >>> ind1 + array([ 0, 9, 32, 53, 74, 95, 116, 137]) + >>> tp[ind1] + array([-0.00743352, 1.08753972, -1.07206545, 1.09550837, -1.07940458, + 1.07849396, -1.0995006 , 1.08094452]) + ''' + + +def test_rfcfilter(): + ''' + # 1. Filtered signal y is the turning points of x. + >>> x = sea() + >>> y = rfcfilter(x[:,1], h=0, method=1) + >>> y[0:5] + array([-1.2004945 , 0.83950546, -0.09049454, -0.02049454, -0.09049454]) + + # 2. This removes all rainflow cycles with range less than 0.5. + >>> y1 = rfcfilter(x[:,1], h=0.5) + >>> y1[0:5] + array([-1.2004945 , 0.83950546, -0.43049454, 0.34950546, -0.51049454]) + + >>> t = linspace(0,7*pi,250) + >>> x = sin(t)+0.1*sin(50*t) + >>> ind = findextrema(x) + >>> ind + array([ 1, 3, 4, 6, 7, 9, 11, 13, 14, 16, 18, 19, 21, + 23, 25, 26, 28, 29, 31, 33, 35, 36, 38, 39, 41, 43, + 45, 46, 48, 50, 51, 53, 55, 56, 58, 60, 61, 63, 65, + 67, 68, 70, 71, 73, 75, 77, 78, 80, 81, 83, 85, 87, + 88, 90, 92, 93, 95, 97, 99, 100, 102, 103, 105, 107, 109, + 110, 112, 113, 115, 117, 119, 120, 122, 124, 125, 127, 129, 131, + 132, 134, 135, 137, 139, 141, 142, 144, 145, 147, 149, 151, 152, + 154, 156, 157, 159, 161, 162, 164, 166, 167, 169, 171, 173, 174, + 176, 177, 179, 181, 183, 184, 186, 187, 189, 191, 193, 194, 196, + 198, 199, 201, 203, 205, 206, 208, 209, 211, 213, 215, 216, 218, + 219, 221, 223, 225, 226, 228, 230, 231, 233, 235, 237, 238, 240, + 241, 243, 245, 247, 248]) + >>> ti, tp = t[ind], x[ind] + >>> tp03 = rfcfilter(tp,0.3) + >>> tp03 + array([-0.00743352, 1.08753972, -1.07206545, 1.09550837, -1.07940458, + 1.07849396, -1.0995006 , 1.08094452, 0.11983423]) + ''' + + +def test_findtp(): + ''' + >>> import numpy as np + >>> x = sea() + >>> x1 = x[0:200,:] + >>> itp = findtp(x1[:,1],0,'Mw') + >>> itph = findtp(x1[:,1],0.3,'Mw') + >>> itp + array([ 11, 21, 22, 24, 26, 28, 31, 39, 43, 45, 47, 51, 56, + 64, 70, 78, 82, 84, 89, 94, 101, 108, 119, 131, 141, 148, + 149, 150, 159, 173, 184, 190, 199]) + >>> itph + array([ 11, 28, 31, 39, 47, 51, 56, 64, 70, 78, 89, 94, 101, + 108, 119, 131, 141, 148, 159, 173, 184, 190, 199]) + ''' + + +def test_findtc(): + ''' + >>> x = sea() + >>> x1 = x[0:200,:] + >>> itc, iv = findtc(x1[:,1],0,'dw') + >>> itc + array([ 28, 31, 39, 56, 64, 69, 78, 82, 83, 89, 94, 101, 108, + 119, 131, 140, 148, 159, 173, 184]) + >>> iv + array([ 19, 29, 34, 53, 60, 67, 76, 81, 82, 84, 90, 99, 103, + 112, 127, 137, 143, 154, 166, 180, 185]) + ''' + + +def test_findoutliers(): + ''' + >>> xx = sea() + >>> dt = diff(xx[:2,0]) + >>> dcrit = 5*dt + >>> ddcrit = 9.81/2*dt*dt + >>> zcrit = 0 + >>> [inds, indg] = findoutliers(xx[:,1],zcrit,dcrit,ddcrit,verbose=True) + Found 0 spurious positive jumps of Dx + Found 0 spurious negative jumps of Dx + Found 37 spurious positive jumps of D^2x + Found 200 spurious negative jumps of D^2x + Found 244 consecutive equal values + Found the total of 1152 spurious points + >>> inds + array([ 6, 7, 8, ..., 9509, 9510, 9511]) + >>> indg + array([ 0, 1, 2, ..., 9521, 9522, 9523]) + ''' + + +def test_common_shape(): + ''' + >>> import numpy as np + >>> A = np.ones((4,1)) + >>> B = 2 + >>> C = np.ones((1,5))*5 + >>> common_shape(A,B,C) + (4, 5) + >>> common_shape(A,B,C,shape=(3,4,1)) + (3, 4, 5) + >>> A = np.ones((4,1)) + >>> B = 2 + >>> C = np.ones((1,5))*5 + >>> common_shape(A,B,C) + (4, 5) + >>> common_shape(A,B,C,shape=(3,4,1)) + (3, 4, 5) + ''' + + +def test_argsreduce(): + ''' + >>> import numpy as np + >>> rand = np.random.random_sample + >>> A = linspace(0,19,20).reshape((4,5)) + >>> B = 2 + >>> C = range(5) + >>> cond = np.ones(A.shape) + >>> [A1,B1,C1] = argsreduce(cond,A,B,C) + >>> B1.shape + (20,) + >>> cond[2,:] = 0 + >>> [A2,B2,C2] = argsreduce(cond,A,B,C) + >>> B2.shape + (15,) + >>> A2;B2;C2 + array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 15., + 16., 17., 18., 19.]) + array([2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]) + array([0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4]) + ''' + + +def test_stirlerr(): + ''' + >>> stirlerr(range(5)) + array([ inf, 0.08106147, 0.0413407 , 0.02767793, 0.02079067]) + ''' + + +def test_getshipchar(): + ''' + >>> sc = getshipchar(10,'service_speed') + >>> names = ['beam', 'beamSTD', 'draught', + ... 'draughtSTD', 'length', 'lengthSTD', + ... 'max_deadweight', 'max_deadweightSTD', 'propeller_diameter', + ... 'propeller_diameterSTD', 'service_speed', 'service_speedSTD'] + >>> for name in names: print( '%s : %g' % (name, sc[name])) + beam : 29 + beamSTD : 2.9 + draught : 9.6 + draughtSTD : 2.112 + length : 216 + lengthSTD : 2.01131 + max_deadweight : 30969 + max_deadweightSTD : 3096.9 + propeller_diameter : 6.76117 + propeller_diameterSTD : 0.20267 + service_speed : 10 + service_speedSTD : 0 + ''' + + +def test_betaloge(): + ''' + >>> betaloge(3, arange(4)) + array([ inf, -1.09861229, -2.48490665, -3.40119738]) + ''' + + +def test_gravity(): + ''' + >>> phi = linspace(0,45,5) + >>> gravity(phi) + array([ 9.78049 , 9.78245014, 9.78803583, 9.79640552, 9.80629387]) + ''' + + +def test_nextpow2(): + ''' + >>> nextpow2(10) + 4 + >>> nextpow2(np.arange(5)) + 3 + ''' + + +def test_discretize(): + ''' + >>> x, y = discretize(np.cos,0,np.pi) + >>> x; y + array([ 0. , 0.19634954, 0.39269908, 0.58904862, 0.78539816, + 0.9817477 , 1.17809725, 1.37444679, 1.57079633, 1.76714587, + 1.96349541, 2.15984495, 2.35619449, 2.55254403, 2.74889357, + 2.94524311, 3.14159265]) + array([ 1.00000000e+00, 9.80785280e-01, 9.23879533e-01, + 8.31469612e-01, 7.07106781e-01, 5.55570233e-01, + 3.82683432e-01, 1.95090322e-01, 6.12323400e-17, + -1.95090322e-01, -3.82683432e-01, -5.55570233e-01, + -7.07106781e-01, -8.31469612e-01, -9.23879533e-01, + -9.80785280e-01, -1.00000000e+00]) + ''' + + +def test_discretize_adaptive(): + ''' + >>> x, y = discretize(np.cos,0,np.pi, method='adaptive') + >>> x; y + array([ 0. , 0.19634954, 0.39269908, 0.58904862, 0.78539816, + 0.9817477 , 1.17809725, 1.37444679, 1.57079633, 1.76714587, + 1.96349541, 2.15984495, 2.35619449, 2.55254403, 2.74889357, + 2.94524311, 3.14159265]) + array([ 1.00000000e+00, 9.80785280e-01, 9.23879533e-01, + 8.31469612e-01, 7.07106781e-01, 5.55570233e-01, + 3.82683432e-01, 1.95090322e-01, 6.12323400e-17, + -1.95090322e-01, -3.82683432e-01, -5.55570233e-01, + -7.07106781e-01, -8.31469612e-01, -9.23879533e-01, + -9.80785280e-01, -1.00000000e+00]) + ''' + + +def test_pol2cart_n_cart2pol(): + ''' + >>> r = 5 + >>> t = linspace(0,pi,20) + >>> x, y = polar2cart(t,r) + >>> x; y + array([ 5. , 4.93180652, 4.72908621, 4.39736876, 3.94570255, + 3.38640786, 2.73474079, 2.00847712, 1.22742744, 0.41289673, + -0.41289673, -1.22742744, -2.00847712, -2.73474079, -3.38640786, + -3.94570255, -4.39736876, -4.72908621, -4.93180652, -5. ]) + array([ 0.00000000e+00, 8.22972951e-01, 1.62349735e+00, + 2.37973697e+00, 3.07106356e+00, 3.67861955e+00, + 4.18583239e+00, 4.57886663e+00, 4.84700133e+00, + 4.98292247e+00, 4.98292247e+00, 4.84700133e+00, + 4.57886663e+00, 4.18583239e+00, 3.67861955e+00, + 3.07106356e+00, 2.37973697e+00, 1.62349735e+00, + 8.22972951e-01, 6.12323400e-16]) + >>> ti, ri = cart2polar(x,y) + >>> ti;ri + array([ 0. , 0.16534698, 0.33069396, 0.49604095, 0.66138793, + 0.82673491, 0.99208189, 1.15742887, 1.32277585, 1.48812284, + 1.65346982, 1.8188168 , 1.98416378, 2.14951076, 2.31485774, + 2.48020473, 2.64555171, 2.81089869, 2.97624567, 3.14159265]) + array([ 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., + 5., 5., 5., 5., 5., 5., 5.]) + ''' + + +def test_meshgrid(): + ''' + >>> x = np.linspace(0,1,3) # coordinates along x axis + >>> y = np.linspace(0,1,2) # coordinates along y axis + >>> xv, yv = meshgrid(x,y) # extend x and y for a 2D xy grid + >>> xv + array([[ 0. , 0.5, 1. ], + [ 0. , 0.5, 1. ]]) + >>> yv + array([[ 0., 0., 0.], + [ 1., 1., 1.]]) + >>> xv, yv = meshgrid(x,y, sparse=True) # make sparse output arrays + >>> xv + array([[ 0. , 0.5, 1. ]]) + >>> yv + array([[ 0.], + [ 1.]]) + + >>> meshgrid(x,y,sparse=True,indexing='ij') # change to matrix indexing + [array([[ 0. ], + [ 0.5], + [ 1. ]]), array([[ 0., 1.]])] + >>> meshgrid(x,y,indexing='ij') + [array([[ 0. , 0. ], + [ 0.5, 0.5], + [ 1. , 1. ]]), array([[ 0., 1.], + [ 0., 1.], + [ 0., 1.]])] + + >>> meshgrid(0,1,5) # just a 3D point + [array([[[0]]]), array([[[1]]]), array([[[5]]])] + >>> map(np.squeeze,meshgrid(0,1,5)) # just a 3D point + [array(0), array(1), array(5)] + >>> meshgrid(3) + array([3]) + >>> meshgrid(y) # 1D grid y is just returned + array([ 0., 1.]) + + `meshgrid` is very useful to evaluate functions on a grid. + + >>> x = np.arange(-5, 5, 0.1) + >>> y = np.arange(-5, 5, 0.1) + >>> xx, yy = meshgrid(x, y, sparse=True) + >>> z = np.sin(xx**2+yy**2)/(xx**2+yy**2) + ''' + + +def test_tranproc(): + ''' + >>> import wafo.transform.models as wtm + >>> tr = wtm.TrHermite() + >>> x = linspace(-5,5,501) + >>> g = tr(x) + >>> y0, y1 = tranproc(x, g, range(5), ones(5)) + >>> y0;y1 + array([ 0.02659612, 1.00115284, 1.92872532, 2.81453257, 3.66292878]) + array([ 1.00005295, 0.9501118 , 0.90589954, 0.86643821, 0.83096482]) + ''' +if __name__ == '__main__': + import doctest + doctest.testmod()