Added wafo.stats + small fixes

pywafo/src/wafo/__init__.py

@@ -7,3 +7,4 @@ import wafo.spectrum
 import wafo.transform
 import wafo.definitions
 import wafo.polynomial
+import wafo.stats

pywafo/src/wafo/misc.py

@@ -101,14 +101,14 @@ class Bunch(object):
     ''' Implement keyword argument initialization of class
 
     '''
-    def __init__(self, ** kwargs):
+    def __init__(self, **kwargs):
         self.__dict__.update(kwargs)
     def keys(self):
         return self.__dict__.keys()
     def update(self, ** kwargs):
         self.__dict__.update(kwargs)
 
-def printf(format, * args):
+def printf(format, *args):
     sys.stdout.write(format % args)
 
 
@@ -135,7 +135,7 @@ def sub_dict_select(somedict, somekeys):
     return dict((k, somedict[k]) for k in somekeys if k in somedict)
 
 
-def parse_kwargs(options, ** kwargs):
+def parse_kwargs(options, **kwargs):
     ''' Update options dict from keyword arguments if the keyword exists in options
 
     Example
@@ -154,11 +154,11 @@ def parse_kwargs(options, ** kwargs):
         options.update(newopts)
     return options
 
-def testfun(*args, ** kwargs):
+def testfun(*args, **kwargs):
     opts = dict(opt1=1, opt2=2)
     if len(args) == 1 and len(kwargs) == 0 and type(args[0]) is str and args[0].startswith('default'):
         return opts
-    opts = parse_kwargs(opts, ** kwargs)
+    opts = parse_kwargs(opts, **kwargs)
     return opts
 
 def detrendma(x, L):
@@ -261,6 +261,9 @@ def ecross(t, f, ind, v):
     --------
     findcross
     '''
+    # Tested on: Python 2.5
+    # revised pab Feb2004
+    # By pab 18.06.2001
     return t[ind] + (v - f[ind]) * (t[ind + 1] - t[ind]) / (f[ind + 1] - f[ind])
 
 def _findcross(xn):

pywafo/src/wafo/objects.py

@@ -47,6 +47,7 @@ __all__ = ['TimeSeries', 'LevelCrossings', 'CyclePairs', 'TurningPoints',
 
 
 
+
 class LevelCrossings(WafoData):
     '''
     Container class for Level crossing data objects in WAFO
@@ -105,6 +106,7 @@ class LevelCrossings(WafoData):
         Example
         -------
         >>> import wafo.spectrum.models as sm
+        >>> from wafo.objects import mat2timeseries
         >>> Sj = sm.Jonswap(Hm0=7)
         >>> S = Sj.tospecdata()   #Make spectrum object from numerical values
         >>> alpha = S.characteristic('alpha')[0]
@@ -243,7 +245,7 @@ class LevelCrossings(WafoData):
             Smaller values gives smoother functions.
         param : 
             vector which defines the region of variation of the data X.
-                     (default [-5 5 513]).    
+                     (default [-5, 5, 513]).    
         monitor : bool
             if true monitor development of estimation
         linextrap : bool
@@ -259,9 +261,9 @@ class LevelCrossings(WafoData):
         ntr :  scalar integer
             Maximum length of empirical crossing intensity. The empirical 
             crossing intensity is interpolated linearly  before smoothing if the
-             length exceeds ntr. A reasonable NTR will significantly speed up the
+            length exceeds ntr. A reasonable NTR (eg. 1000) will significantly 
-            estimation for long time series without loosing any accuracy. 
+            speed up the estimation for long time series without loosing any accuracy. 
-            NTR should be chosen greater than PARAM(3). (default 1000)
+            NTR should be chosen greater than PARAM(3). (default inf)
             
         Returns
         -------
@@ -284,11 +286,12 @@ class LevelCrossings(WafoData):
         -------
         >>> import wafo.spectrum.models as sm
         >>> import wafo.transform.models as tm
+        >>> from wafo.objects import mat2timeseries
         >>> Hs = 7.0
         >>> Sj = sm.Jonswap(Hm0=Hs)
         >>> S = Sj.tospecdata()   #Make spectrum object from numerical values
-        >>> S.tr = tm.TrOchi(mean=0, skew=sk, kurt=0, sigma=Hs/4, ysigma=Hs/4)
+        >>> S.tr = tm.TrOchi(mean=0, skew=0.16, kurt=0, sigma=Hs/4, ysigma=Hs/4)
-        >>> xs = S.sim(ns=2**13)
+        >>> xs = S.sim(ns=2**16)
         >>> ts = mat2timeseries(xs)
         >>> tp = ts.turning_points()
         >>> mm = tp.cycle_pairs()
@@ -296,6 +299,17 @@ class LevelCrossings(WafoData):
         >>> g0, gemp = lc.trdata(monitor=True) # Monitor the development
         >>> g1, gemp = lc.trdata(gvar=0.5 ) # Equal weight on all points
         >>> g2, gemp = lc.trdata(gvar=[3.5, 0.5, 3.5])  # Less weight on the ends
+        >>> S.tr.dist2gauss()
+        5.9322684525265501
+        >>> np.round(gemp.dist2gauss())
+        6.0
+        >>> np.round(g0.dist2gauss())
+        4.0
+        >>> np.round(g1.dist2gauss())
+        4.0
+        >>> np.round(g2.dist2gauss())
+        4.0
+        
          hold on, trplot(g1,g)                          # Check the fit
          trplot(g2)
         
@@ -324,7 +338,7 @@ class LevelCrossings(WafoData):
             sigma = self.stdev
         
         opt = DotDict(chkder=True, plotflag=False, csm=0.9, gsm=.05,
-            param=(-5, 5, 513), delay=2, lin_extrap=True, ntr=1000, ne=7, cvar=1, gvar=1)
+            param=(-5, 5, 513), delay=2, lin_extrap=True, ntr=inf, ne=7, cvar=1, gvar=1)
         # If just 'defaults' passed in, return the default options in g
         
         opt.update(options)
@@ -341,13 +355,13 @@ class LevelCrossings(WafoData):
         else:
             Ner = 0
             lc1, lc2 = self.args, self.data        
-        ng = len(opt.gvar)
+        ng = len(atleast_1d(opt.gvar))
         if ng == 1:
             gvar = opt.gvar * ones(ncr)
         else:
             gvar = interp1d(linspace(0, 1, ng) , opt.gvar, kind='linear')( linspace(0, 1, ncr))  
         
-        ng = len(opt.cvar)
+        ng = len(atleast_1d(opt.cvar))
         if ng == 1:
             cvar = opt.cvar * ones(ncr)
         else:
@@ -369,7 +383,7 @@ class LevelCrossings(WafoData):
             cor2 = 0
         
     
-        lc22 = cumtrapz(lc2, lc1) + cor1
+        lc22 = hstack((0,cumtrapz(lc2, lc1) + cor1))
         lc22 = (lc22 + 0.5) / (lc22[-1] + cor2 + 1)
         lc11 = (lc1 - mean) / sigma
         
@@ -418,7 +432,7 @@ class LevelCrossings(WafoData):
         if opt.plotflag > 0:
             g.plot()
             g2.plot()
-            
+        g2.setplotter('step')
         return g, g2
 
 class CyclePairs(WafoData):
@@ -679,7 +693,7 @@ class TimeSeries(WafoData):
     --------
     >>> import wafo.data
     >>> x = wafo.data.sea()
-    >>> ts = mat2timeseries(x)
+    >>> ts = wafo.objects.mat2timeseries(x)
     >>> rf = ts.tocovdata(lag=150)
     >>> h = rf.plot()
 
@@ -688,7 +702,7 @@ class TimeSeries(WafoData):
         super(TimeSeries, self).__init__(*args, **kwds)
         self.name = 'WAFO TimeSeries Object'
         self.sensortypes = ['n', ]
-        self.position = zeros(3)
+        self.position = [zeros(3), ]
         somekeys = ['sensortypes', 'position']
         self.__dict__.update(sub_dict_select(kwds, somekeys))
 

pywafo/src/wafo/spectrum/core.py

@@ -1,4 +1,5 @@
 from __future__ import division
+from scipy.misc.ppimport import ppimport
 import warnings
 import numpy as np
 from numpy import (pi, inf, meshgrid, zeros, ones, where, nonzero, #@UnresolvedImport
@@ -32,6 +33,7 @@ from wafo.plotbackend import plotbackend
 
 
 # Trick to avoid error due to circular import
+#_WAFOCOV = ppimport('wafo.covariance')
 _WAFOCOV = JITImport('wafo.covariance')
 
 
@@ -366,7 +368,7 @@ class SpecData1D(WafoData):
             txt = '''Spectrum does not start at zero frequency/wave number.
             Correct it with resample, for example.'''
             raise ValueError(txt)
-        d_w = abs(diff(freq, n_f=2, axis=0))
+        d_w = abs(diff(freq, n=2, axis=0))
         if any(d_w > 1.0e-8):
             txt = '''Not equidistant frequencies/wave numbers in spectrum. 
             Correct it with resample, for example.'''
@@ -390,7 +392,7 @@ class SpecData1D(WafoData):
         if self.freqtype in 'k':
             lagtype = 'x'
         else:
-            lagtype = 'time'
+            lagtype = 't'
 
         d_t = spec.sampling_period()
         #normalize spec so that sum(specn)/(n_f-1)=acf(0)=var(X)
@@ -414,7 +416,7 @@ class SpecData1D(WafoData):
 
         if nr > 0:
             w = r_[w , zeros(nfft - 2 * n_f + 2) , -w[n_f - 1:0:-1] ]
-            fieldname = 'R' + lagtype * nr 
+            fieldname = 'R' + lagtype[0] * nr 
             for i in range(1, nr + 1):
                 rper = -1j * w * rper
                 d_acf = fft(rper, nfft).real / (2 * n_f - 2)
@@ -1481,6 +1483,108 @@ class SpecData1D(WafoData):
 ##            skew = sum((6*C2+8*E2).*E)/sa^3   % skewness
 ##            kurt = 3+48*sum((C2+E2).*E2)/sa^4 % kurtosis
         return output
+    def testgaussian(ns,test0=None, cases=100, method='nonlinear',**opt):
+        '''
+        TESTGAUSSIAN Test if a stochastic process is Gaussian.
+        
+         CALL:  test1 = testgaussian(S,[ns,Ns],test0,def,options);
+        
+         test1,
+            test0 = simulated and observed value of e(g)=int (g(u)-u)^2 du,
+                    respectively, where int limits is given by OPTIONS.PARAM. 
+                 
+              S   = spectral density structure 
+               ns = # of points simulated
+               cases = # of independent simulations (default  100)
+        
+           def    = 'nonlinear' : transform based on smoothed crossing intensity (default)
+                    'mnonlinear': transform based on smoothed marginal distribution
+          options = options structure defining how the estimation of the
+                    transformation is done. (default troptset('dat2tr'))
+        
+         TESTGAUSSIAN simulates  e(g(u)-u) = int (g(u)-u)^2 du  for Gaussian processes 
+         given the spectral density, S. The result is plotted if test0 is given.
+         This is useful for testing if the process X(t) is Gaussian.
+         If 95% of TEST1 is less than TEST0 then X(t) is not Gaussian at a 5% level.
+        
+         Example:
+         Hm0 = 7;
+         S0 = jonswap([],Hm0); g=ochitr([],[Hm0/4]); S=S0;
+         S.tr=g;S.tr(:,2)=g(:,2)*Hm0/4;
+         xs = spec2sdat(S,2^13);
+         [g0 t0] = dat2tr(xs);
+         t1 = testgaussian(S0,[2^13 50],t0); 
+        
+         See also  cov2sdat, dat2tr, troptset
+        '''
+#        Tested on: Matlab 5.3, 5.2, 5.1
+#         History:
+#         revised pab
+#         -changed name from mctrtest to testgaussianity
+#         revised pab jan2005
+#         changed order of input so that chapter1.m works
+#         revised pab Feb2004
+#         -changed h1line  
+#         revised pab 29.12.2000
+#         - added options and def to input due to changed calling syntax for dat2tr.
+#         - updated the help header
+#         revised by pab 12.11.99
+#            fixed a bug, string input to dat2tr 'nonlinear' 
+#         revised by pab 12.10.99
+#            updated help header 
+#         revised by pab 11.08.99
+#         changed name from mctest to mctrtest
+#         by pab 11.11.98
+        
+        maxsize = 200000 # must divide the computations due to limited memory
+#        if nargin<5||isempty(opt):
+#            opt = troptset('dat2tr');
+#        
+#        opt = troptset(opt,'multip',1)
+        
+        if test0 is None:
+            plotflag=0
+        else: 
+            plotflag=1
+        
+        if cases>50:
+            print('  ... be patient this may take a while')
+        
+        test1 = []
+        rep = floor(ns*cases/maxsize)+1
+        
+        Nstep = floor(cases/rep);
+        
+        acf = self.tocovdata()
+        #R = spec2cov(S);
+        
+        for ix in  range(rep):
+            xs = acf.sim(ns=ns, cases=Nstep)
+            #xs = cov2sdat(R,[ns Nstep]);
+            [g, tmp] = dat2tr(xs,method, **opt);
+            #test1 = [test1; tmp(:)]
+            print('finished %d of %d ' % (ix,rep) )
+        
+        if rep>1:
+            xs = acf.sim(ns=ns, cases=rem(cases,rep))
+            [g, tmp] = dat2tr(xs,method,**opt);
+            #test1 = [test1; tmp(:)];
+        
+        
+#        if (nargout>0 || plotflag==0),
+#          test2=test1;
+#        end
+#        
+#        
+#        if plotflag 
+#          plot(test1,'o'),hold on
+#          if 1 
+#            plot([1 cases],test0*[1 1],'--'),
+#          end
+#          hold off
+#          ylabel('e(g(u)-u)')
+#          xlabel('Simulation number')
+#        end
         
     def moment(self, nr=2, even=True, j=0):
         ''' Calculates spectral moments from spectrum

pywafo/src/wafo/stats/__init__.py

@@ -0,0 +1,13 @@
+"""
+Statistics package in WAFO Toolbox.
+
+    Readme               - Readme file for module STATS in WAFO Toolbox
+
+"""
+from scipy.stats import *
+from wafo.stats.core import *
+#from wafo.stats.distributions import *
+#from wafo.spectrum.core import SpecData1D
+#import wafo.spectrum.models
+#import wafo.spectrum.dispersion_relation
+#from wafo.data_structures import SpecData1D

pywafo/src/wafo/stats/core.py

@@ -0,0 +1,86 @@
+from __future__ import division
+from wafo.wafodata import WafoData
+import numpy as np
+from numpy import inf
+from numpy.lib.shape_base import atleast_1d
+from numpy.ma.core import arange, floor
+__all__ = ['edf']
+
+def edf(x, method=2):
+    ''' 
+    Returns EDF Empirical Distribution Function.
+    
+    Parameters
+    ----------
+    x : array-like
+        data vector
+    method : integer scalar 
+        1. Interpolation so that F(X_(k)) == (k-0.5)/n.
+        2. Interpolation so that F(X_(k)) == k/(n+1).    (default)
+        3. The empirical distribution. F(X_(k)) = k/n
+     
+    Example
+    -------
+    >>> import wafo.stats as ws
+    >>> x = np.linspace(0,6,200)
+    >>> R = ws.rayleigh.rvs(scale=2,size=100)
+    >>> F = ws.edf(R)
+    >>> F.plot()
+      
+     See also edf, pdfplot, cumtrapz
+    '''
+        
+
+    z = atleast_1d(x)       
+    z.sort()
+    
+    N = len(z)
+    if method == 1:
+        Fz1 = arange(0.5, N) / N
+    elif method == 3:
+        Fz1 = arange(1, N + 1) / N
+    else:
+        Fz1 = arange(1, N + 1) / (N + 1)
+     
+    F = WafoData(Fz1, z, xlab='x', ylab='F(x)')
+    F.setplotter('step')
+    return F
+
+def edfcnd(x, c=None, method=2):
+    ''' 
+    EDFCND Empirical Distribution Function CoNDitioned that X>=c.
+    
+    Parameters
+    ----------
+    x : array-like
+        data vector
+    method : integer scalar 
+        1. Interpolation so that F(X_(k)) == (k-0.5)/n.
+        2. Interpolation so that F(X_(k)) == k/(n+1).    (default)
+        3. The empirical distribution. F(X_(k)) = k/n
+     
+    Example
+    -------
+    >>> import wafo.stats as ws
+    >>> x = np.linspace(0,6,200)
+    >>> R = ws.rayleigh.rvs(scale=2,size=100)
+    >>> Fc = ws.edfcd(R, 1)
+    >>> Fc.plot()
+    >>> F = ws.edf(R)
+    >>> F.plot()
+    
+     See also edf, pdfplot, cumtrapz
+    '''
+    z = atleast_1d(x)
+    if c is None:
+        c = floor(min(z.min(), 0))
+    
+    try:
+        F = edf(z[c <= z], method=method)
+    except:
+        ValueError('No data points above c=%d' % int(c)) 
+    
+    if - inf < c:
+        F.labels.ylab = 'F(x| X>=%g)' % c
+    
+    return F

pywafo/src/wafo/stats/misc.py

@@ -0,0 +1,97 @@
+from numpy import asarray, ndarray, reshape, repeat, nan, product
+
+def valarray(shape,value=nan,typecode=None):
+    """Return an array of all value.
+    """
+    out = reshape(repeat([value],product(shape,axis=0),axis=0),shape)
+    if typecode is not None:
+        out = out.astype(typecode)
+    if not isinstance(out, ndarray):
+        out = asarray(out)
+    return out
+
+class rv_frozen(object):
+    ''' Frozen continous or discrete 1D Random Variable object (RV)
+
+    Methods
+    -------
+    RV.rvs(size=1)
+        - random variates
+
+    RV.pdf(x)
+        - probability density function (continous case)
+
+    RV.pmf(x)
+        - probability mass function (discrete case)
+
+    RV.cdf(x)
+        - cumulative density function
+
+    RV.sf(x)
+        - survival function (1-cdf --- sometimes more accurate)
+
+    RV.ppf(q)
+        - percent point function (inverse of cdf --- percentiles)
+
+    RV.isf(q)
+        - inverse survival function (inverse of sf)
+
+    RV.stats(moments='mv')
+        - mean('m'), variance('v'), skew('s'), and/or kurtosis('k')
+
+    RV.entropy()
+        - (differential) entropy of the RV.
+
+    Parameters
+    ----------
+    x : array-like
+        quantiles
+    q : array-like
+        lower or upper tail probability
+    size : int or tuple of ints, optional, keyword
+        shape of random variates
+    moments : string, optional, keyword
+        one or more of 'm' mean, 'v' variance, 's' skewness, 'k' kurtosis
+    '''
+    def __init__(self, dist, *args, **kwds):
+        self.dist = dist
+        loc0, scale0 = map(kwds.get, ['loc', 'scale'])
+        if isinstance(dist,rv_continuous):
+            args, loc0, scale0 = dist.fix_loc_scale(args, loc0, scale0)
+            self.par = args + (loc0, scale0)
+        else: # rv_discrete
+            args, loc0 = dist.fix_loc(args, loc0)
+            self.par = args + (loc0,)
+
+
+    def pdf(self,x):
+        ''' Probability density function at x of the given RV.'''
+        return self.dist.pdf(x,*self.par)
+    def cdf(self,x):
+        '''Cumulative distribution function at x of the given RV.'''
+        return self.dist.cdf(x,*self.par)
+    def ppf(self,q):
+        '''Percent point function (inverse of cdf) at q of the given RV.'''
+        return self.dist.ppf(q,*self.par)
+    def isf(self,q):
+        '''Inverse survival function at q of the given RV.'''
+        return self.dist.isf(q,*self.par)
+    def rvs(self, size=None):
+        '''Random variates of given type.'''
+        kwds = dict(size=size)
+        return self.dist.rvs(*self.par,**kwds)
+    def sf(self,x):
+        '''Survival function (1-cdf) at x of the given RV.'''
+        return self.dist.sf(x,*self.par)
+    def stats(self,moments='mv'):
+        ''' Some statistics of the given RV'''
+        kwds = dict(moments=moments)
+        return self.dist.stats(*self.par,**kwds)
+    def moment(self,n):
+        par1 = self.par[:self.dist.numargs]
+        return self.dist.moment(n,*par1)
+    def entropy(self):
+        return self.dist.entropy(*self.par)
+    def pmf(self,k):
+        '''Probability mass function at k of the given RV'''
+        return self.dist.pmf(k,*self.par)

pywafo/src/wafo/stats/plotbackend.py

@@ -0,0 +1,9 @@
+# Set this to eg. pylab to be able to plot
+import numpy
+try:
+    from matplotlib import pyplot as plotbackend
+    #from matplotlib import pyplot
+    numpy.disp('Scipy.stats: plotbackend is set to matplotlib.pyplot')
+except:
+    numpy.disp('Scipy.stats: Unable to load matplotlib.pyplot as plotbackend')
+    plotbackend = None

pywafo/src/wafo/wafodata.py

@@ -62,7 +62,7 @@ class WafoData(object):
     specdata,
     covdata
     '''
-    def __init__(self, data=None, args=None,**kwds):
+    def __init__(self, data=None, args=None,*args2,**kwds):
         self.data = data
         self.args = args
         self.date = now()