diff --git a/pywafo/src/wafo/objects.py b/pywafo/src/wafo/objects.py index f303e1e..6a5a36b 100644 --- a/pywafo/src/wafo/objects.py +++ b/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,12 +2403,8 @@ 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): ''' Class for computing transfer functions based on linear wave theory @@ -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 diff --git a/pywafo/src/wafo/stats/distributions.py b/pywafo/src/wafo/stats/distributions.py index 8d5daee..5035251 100644 --- a/pywafo/src/wafo/stats/distributions.py +++ b/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() diff --git a/pywafo/src/wafo/transform/core.py b/pywafo/src/wafo/transform/core.py index abebfa6..47d5cfc 100644 --- a/pywafo/src/wafo/transform/core.py +++ b/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) diff --git a/pywafo/src/wafo/transform/models.py b/pywafo/src/wafo/transform/models.py index d59afd5..4134645 100644 --- a/pywafo/src/wafo/transform/models.py +++ b/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)