Small bugfixes and cosmetic fixes.

pywafo/src/wafo/objects.py

@@ -167,12 +167,12 @@ class LevelCrossings(WafoData):
         >>> lc_exp = lc.extrapolate(-2*s, 2*s, dist='expon')
         >>> lc_ray = lc.extrapolate(-2*s, 2*s, dist='rayleigh')
         
+        lc.plot()
         lc_gpd.plot()
         lc_exp.plot()
         lc_ray.plot()
         
         
-        
         See also  
         --------
         cmat2extralc, rfmextrapolate, lc2rfmextreme, extralc, fitgenpar
@@ -246,7 +246,7 @@ class LevelCrossings(WafoData):
         if dist.startswith('gen'): 
             genpareto = distributions.genpareto
             phat = genpareto.fit2(x, floc=0, method=method)
-            F = phat.cdf(xF)
+            SF = phat.sf(xF)
             
             covar = phat.par_cov[::2, ::2]
             # Calculate 90 # confidence region, an ellipse, for (k,s)
@@ -274,28 +274,28 @@ class LevelCrossings(WafoData):
             for i in range(Nc):
                 k = c1[0, i]
                 s = c1[1, i]
-                F2 = genpareto.cdf(xF, k, scale=s)
-                lcEstCu[:, i] = (lcu + dXX) * (1 - F2)
-                lcEstCl[:, i] = (lcu - dXX) * (1 - F2)
+                SF2 = genpareto.sf(xF, k, scale=s)
+                lcEstCu[:, i] = (lcu + dXX) * (SF2)
+                lcEstCl[:, i] = (lcu - dXX) * (SF2)
             #end
             
-            lcEst = np.vstack((xF + u, lcu * (1 - F),
+            lcEst = np.vstack((xF + u, lcu * (SF),
                                  lcEstCl.min(axis=1), lcEstCu.max(axis=1))).T            
         elif dist.startswith('exp'):
             expon = distributions.expon
             phat = expon.fit2(x, floc=0, method=method)           
-            F = phat.cdf(xF) 
-            lcEst =  np.vstack((xF + u, lcu * (1 - F))).T
+            SF = phat.sf(xF) 
+            lcEst =  np.vstack((xF + u, lcu * (SF))).T
             
         elif dist.startswith('ray') or dist.startswith('trun'):
             phat = distributions.truncrayleigh.fit2(x, floc=0, method=method)
-            F = phat.cdf(xF)
+            SF = phat.sf(xF)
             if False:
                 n = len(x)
                 Sx = sum((x + offset) ** 2 - offset ** 2)
                 s = sqrt(Sx / n);  # Shape parameter
                 F = -np.expm1(-((xF + offset) ** 2 - offset ** 2) / s ** 2)
-            lcEst = np.vstack((xF + u, lcu * (1 - F))).T
+            lcEst = np.vstack((xF + u, lcu * (SF))).T
         else: 
             raise ValueError()
           
@@ -2393,7 +2393,7 @@ def sensor_type(*sensorids):
 
     See also 
     --------
-    sensortypeid, tran
+    sensor_typeid, tran
     '''
     valid_names = ('n', 'n_t', 'n_tt', 'n_x', 'n_y', 'n_xx', 'n_yy', 'n_xy',
                   'p', 'u', 'v', 'w', 'u_t', 'v_t', 'w_t', 'x_p', 'y_p', 'z_p',
@@ -2403,11 +2403,7 @@ def sensor_type(*sensorids):
         ids = hstack(ids)
     n = len(valid_names) - 1
     ids = where(((ids < 0) | (n < ids)), n , ids)
-    
-    #try:
     return tuple(valid_names[i] for i in ids)
-    #except:
-    #    raise ValueError('Input must be an integer!')
     
 class TransferFunction(object):
     '''
@@ -2520,6 +2516,7 @@ class TransferFunction(object):
                                X_p=self._x_p, x_p=self._x_p,
                                Y_p=self._y_p, y_p=self._y_p,
                                Z_p=self._z_p, z_p=self._z_p)
+    __call__ = tran 
     def tran(self, w, theta=0, kw=None):
         '''
         Return transfer functions based on linear wave theory
@@ -2578,7 +2575,7 @@ class TransferFunction(object):
     
 #---Private member methods
     def _get_ee_cthxy(self, theta, kw):
-        # convert to from angle in degrees to radians
+        # convert from angle in degrees to radians
         bet = self.bet
         thxr = self.thetax * pi / 180
         thyr = self.thetay * pi / 180

pywafo/src/wafo/stats/distributions.py

@@ -282,7 +282,7 @@ _doc_default_example = \
 """Examples
 --------
 >>> import matplotlib.pyplot as plt
->>> from scipy.stats import %(name)s
+>>> from wafo.stats import %(name)s
 >>> numargs = %(name)s.numargs
 >>> [ %(shapes)s ] = [0.9,] * numargs
 >>> rv = %(name)s(%(shapes)s)
@@ -5151,13 +5151,16 @@ class rayleigh_gen(rv_continuous):
             return x - phat[0] / sqrt(-2.0 * logSF)
         else:
             return x - phat[1] * sqrt(-2.0 * logSF)
-
     def _rvs(self):
         return chi.rvs(2,size=self._size)
     def _pdf(self, r):
         return r*exp(-r*r/2.0)
+    def _logpdf(self, r):
+        return log(r) - r*r/2.0
     def _cdf(self, r):
         return - expm1(-r * r / 2.0)
+    def _sf(self, r):
+        return exp(-r * r / 2.0)
     def _ppf(self, q):
         return sqrt(-2 * log1p(-q))
     def _stats(self):
@@ -5178,26 +5181,45 @@ for x >= 0.
                         )
 
 class truncrayleigh_gen(rv_continuous):
+    def _argcheck(self, c):
+        return (c>=0)
     def link(self, x, logSF, phat, ix):
         rv_continuous.link.__doc__
         c = phat[0]
-        if ix == 1:
-            return x - phat[1] / sqrt(-2.0 * logSF)
-        else:
-            return x - phat[2] * sqrt(-2.0 * logSF)
+        if ix == 2:
+            return x - phat[1] / (sqrt(c*c - 2 * logSF) - c)
+        elif ix == 1:
+            return x - phat[2] * (sqrt(c*c - 2 * logSF) - c)
+        elif ix==0:
+            xn = (x - phat[1])/phat[2]
+            return - 2 * logSF / xn - xn / 2.0
+    def _fitstart(self, data, args=None):
+        if args is None:
+            args = (0.0,)*self.numargs
+        return args + self.fit_loc_scale(data, *args)
     def _pdf(self, r, c):
         rc = r+c
         return rc*exp(-(rc*rc-c*c)/2.0)
+    def _logpdf(self, r, c):
+        rc = r+c
+        return log(rc)-(rc*rc-c*c)/2.0
     def _cdf(self, r, c):
         rc = r+c
         return - expm1(-(rc*rc-c*c)/ 2.0)
+    def _logsf(self, r, c):
+        rc = r+c
+        return -(rc*rc-c*c)/ 2.0
+    def _sf(self, r, c):
+        return exp(self._logsf(r, c))
     def _ppf(self, q, c):
         return sqrt(c*c - 2 * log1p(-q)) - c
     def _stats(self, c):
+        # TODO: correct this it is wrong!
         val = 4-pi
         return np.sqrt(pi/2), val/2, 2*(pi-3)*sqrt(pi)/val**1.5, \
                6*pi/val-16/val**2
     def _entropy(self, c):
+        # TODO: correct this it is wrong!
         return _EULER/2.0 + 1 - 0.5*log(2)
 truncrayleigh = truncrayleigh_gen(a=0.0, name="truncrayleigh",
                         longname="A truncated Rayleigh",
@@ -5207,7 +5229,7 @@ truncrayleigh = truncrayleigh_gen(a=0.0, name="truncrayleigh",
 Truncated Rayleigh distribution
 
 truncrayleigh.cdf(r) = 1 - exp(-((r+c)**2-c**2)/2)
-for x >= 0.
+for x >= 0, c>=0.
 """
                         )
 
@@ -7352,6 +7374,31 @@ Skellam distribution
                       )
                         
 
+def test_truncrayleigh():
+    import matplotlib.pyplot as plt
+    from wafo.stats import truncrayleigh
+    numargs = truncrayleigh.numargs
+    [ c ] = [0.9,] * numargs
+    rv = truncrayleigh(c)
+
+#Display frozen pdf
+
+    x = np.linspace(0, np.minimum(rv.dist.b, 3))
+    h = plt.plot(x, rv.pdf(x))
+
+#Check accuracy of cdf and ppf
+
+    prb = truncrayleigh.cdf(x, c)
+    h = plt.semilogy(np.abs(x - truncrayleigh.ppf(prb, c)) + 1e-20)
+
+#Random number generation
+
+    R = truncrayleigh.rvs(c, size=100)
+
+#Compare ML and MPS method
+    phat = truncrayleigh.fit2(R, method='ml');
+
+ 
 def main():
     import matplotlib
     matplotlib.interactive(True)
@@ -7402,4 +7449,5 @@ def main():
     pht = genpareto.fit(r, 1, par_fix=[0, 0, nan])
     lp = pht.profile()
 if __name__ == '__main__':
-    main()
+    #main()
+    test_truncrayleigh()

pywafo/src/wafo/transform/core.py

@@ -188,6 +188,9 @@ class TrData(WafoData, TrCommon):
             
         self.children = [WafoData((self.args-self.mean)/self.sigma, self.args)]
 
+    def trdata(self):
+        return self
+    
     def _gauss2dat(self, y, *yi):
         return tranproc(self.data, self.args, y, *yi)
     

pywafo/src/wafo/transform/models.py

@@ -44,7 +44,7 @@ _example = '''
     '''
 class TrCommon2(TrCommon):
     __doc__ = TrCommon.__doc__
-    def trdata(self,x=None, xnmin=-5, xnmax=5, n=513):
+    def trdata(self, x=None, xnmin= -5, xnmax=5, n=513):
         """
         Return a discretized transformation model.
 
@@ -66,11 +66,11 @@ class TrCommon2(TrCommon):
         """
         if x is None:
             xn = np.linspace(xnmin, xnmax, n)
-            x = self.sigma*xn+self.mean
+            x = self.sigma * xn + self.mean
         else:
-            xn = (x-self.mean)/self.sigma
+            xn = (x - self.mean) / self.sigma
 
-        yn = (self._dat2gauss(x)-self.ymean)/self.ysigma
+        yn = (self._dat2gauss(x) - self.ymean) / self.ysigma
         
         return TrData(yn, x, mean=self.mean, sigma=self.sigma)