From 9a87ed63b7a3a90e8c08c3445bfaf6a2dad285bf Mon Sep 17 00:00:00 2001 From: pbrod Date: Thu, 9 Jun 2016 03:45:55 +0200 Subject: [PATCH] Small updates --- wafo/spectrum/core.py | 10 +++++++-- wafo/stats/__init__.py | 1 + wafo/stats/_distn_infrastructure.py | 34 +++++++++++++++++++++++++++-- wafo/stats/estimation.py | 20 +++++++++-------- 4 files changed, 52 insertions(+), 13 deletions(-) diff --git a/wafo/spectrum/core.py b/wafo/spectrum/core.py index baee24d..de82919 100644 --- a/wafo/spectrum/core.py +++ b/wafo/spectrum/core.py @@ -1422,8 +1422,14 @@ class SpecData1D(PlotData): The joint density of zero separated Max2min cycles in time (a); in space (b); AcAt in time for nonlinear sea model (c): - Hm0=7;Tp=11 - S = jonswap(4*pi/Tp,[Hm0 Tp]) + >>> from wafo.spectrum import models as sm + >>> w = np.linspace(0,3,100) + >>> Sj = sm.Jonswap() + >>> S = Sj.tospecdata() + >>> f = S.to_t_pdf(pdef='Tc', paramt=(0, 10, 51), speed=7) + >>> Hm0 = 7 + >>> Tp = 11 + >>> S = sm.Jonswap(4*pi/Tp,[Hm0 Tp]) Sk = spec2spec(S,'k1d') L0 = spec2mom(S,1) paramu = [sqrt(L0)*[-4 4] 41] diff --git a/wafo/stats/__init__.py b/wafo/stats/__init__.py index e12a20a..0e5f112 100644 --- a/wafo/stats/__init__.py +++ b/wafo/stats/__init__.py @@ -341,6 +341,7 @@ from .core import * from .distributions import * from . import estimation + # remove vonmises_cython from __all__, I don't know why it is included __all__ = [s for s in dir() if not (s.startswith('_') or s.endswith('cython'))] diff --git a/wafo/stats/_distn_infrastructure.py b/wafo/stats/_distn_infrastructure.py index 8bc802c..bbaa418 100644 --- a/wafo/stats/_distn_infrastructure.py +++ b/wafo/stats/_distn_infrastructure.py @@ -289,7 +289,7 @@ def nlogps(self, theta, x): if not self._argcheck(*args) or scale <= 0: return inf x = asarray((x - loc) / scale) - cond0 = (x <= self.a) | (self.b <= x) + cond0 = (x < self.a) | (self.b < x) Nbad = np.sum(cond0) if Nbad > 0: x = argsreduce(~cond0, x)[0] @@ -328,11 +328,40 @@ def nlogps(self, theta, x): return T +def _penalized_nnlf(self, theta, x): + ''' Return negative loglikelihood function, + i.e., - sum (log pdf(x, theta), axis=0) + where theta are the parameters (including loc and scale) + ''' + try: + loc = theta[-2] + scale = theta[-1] + args = tuple(theta[:-2]) + except IndexError: + raise ValueError("Not enough input arguments.") + if not self._argcheck(*args) or scale <= 0: + return inf + x = asarray((x-loc) / scale) + + loginf = log(_XMAX) + + if np.isneginf(self.a).all() and np.isinf(self.b).all(): + Nbad = 0 + else: + cond0 = (x < self.a) | (self.b < x) + Nbad = sum(cond0) + if Nbad > 0: + x = argsreduce(~cond0, x)[0] + + N = len(x) + return self._nnlf(x, *args) + N*log(scale) + Nbad * 100.0 * loginf + + def _reduce_func(self, args, options): # First of all, convert fshapes params to fnum: eg for stats.beta, # shapes='a, b'. To fix `a`, can specify either `f1` or `fa`. # Convert the latter into the former. - kwds = options.copy() + kwds = options # .copy() if self.shapes: shapes = self.shapes.replace(',', ' ').split() for j, s in enumerate(shapes): @@ -574,6 +603,7 @@ rv_continuous.freeze = freeze rv_continuous.link = link rv_continuous._link = _link rv_continuous.nlogps = nlogps +rv_continuous._penalized_nnlf = _penalized_nnlf rv_continuous._reduce_func = _reduce_func rv_continuous.fit = fit rv_continuous.fit2 = fit2 diff --git a/wafo/stats/estimation.py b/wafo/stats/estimation.py index 6fbd769..2e92c13 100644 --- a/wafo/stats/estimation.py +++ b/wafo/stats/estimation.py @@ -134,7 +134,9 @@ class Profile(object): >>> sf_ci = Lsf.get_bounds(alpha=0.2) ''' - def __init__(self, fit_dist, **kwds): + def __init__(self, fit_dist, method=None, **kwds): + if method is None: + method = fit_dist.method.lower() try: i0 = (1 - np.isfinite(fit_dist.par_fix)).argmax() @@ -153,10 +155,10 @@ class Profile(object): [i0, 100, 0.05, None, None, None, None, None])] self.title = '%g%s CI' % (100 * (1.0 - self.alpha), '%') - if fit_dist.method.startswith('ml'): + if method.startswith('ml'): self.ylabel = self.ylabel + 'likelihood' Lmax = fit_dist.LLmax - elif fit_dist.method.startswith('mps'): + elif method.startswith('mps'): self.ylabel = self.ylabel + ' product spacing' Lmax = fit_dist.LPSmax else: @@ -604,13 +606,13 @@ class FitDistribution(rv_frozen): self.dist = dist numargs = dist.numargs - self.method = method + self.method = method.lower() self.alpha = alpha self.par_fix = par_fix self.search = search self.copydata = copydata - if self.method.lower()[:].startswith('mps'): + if self.method.startswith('mps'): self._fitfun = self._nlogps else: self._fitfun = self._nnlf @@ -797,7 +799,7 @@ class FitDistribution(rv_frozen): return np.inf dist = self.dist x = asarray((x - loc) / scale) - cond0 = (x <= dist.a) | (dist.b <= x) + cond0 = (x < dist.a) | (dist.b < x) Nbad = np.sum(cond0) if Nbad > 0: x = argsreduce(~cond0, x)[0] @@ -816,8 +818,8 @@ class FitDistribution(rv_frozen): if any(tie): # TODO : implement this method for treating ties in data: # Assume measuring error is delta. Then compute - # yL = F(xi-delta,theta) - # yU = F(xi+delta,theta) + # yL = F(xi - delta, theta) + # yU = F(xi + delta, theta) # and replace # logDj = log((yU-yL)/(r-1)) for j = i+1,i+2,...i+r-1 @@ -1020,7 +1022,7 @@ class FitDistribution(rv_frozen): fixstr = 'Fixed: phat[%s] = %s ' % (phatistr, phatvstr) infostr = 'Fit method: %s, Fit p-value: %2.2f %s' % ( - self.method, self.pvalue, fixstr) + self.method.upper(), self.pvalue, fixstr) try: plotbackend.figtext(0.05, 0.01, infostr) except: