diff --git a/wafo/stats/_continuous_distns.py b/wafo/stats/_continuous_distns.py index cdce229..dc7980d 100644 --- a/wafo/stats/_continuous_distns.py +++ b/wafo/stats/_continuous_distns.py @@ -1120,15 +1120,6 @@ class expon_gen(rv_continuous): %(example)s """ - - def _link(self, x, logSF, phat, ix): - if ix == 1: - return - (x - phat[0]) / logSF - elif ix == 0: - return x + phat[1] * logSF - else: - raise IndexError('Index to the fixed parameter is out of bounds') - def _rvs(self): return self._random_state.standard_exponential(self._size) @@ -1593,17 +1584,6 @@ class frechet_r_gen(rv_continuous): %(example)s """ - def _link(self, x, logSF, phat, ix): - if ix == 0: - phati = log(-logSF) / log((x - phat[1]) / phat[2]) - elif ix == 1: - phati = x - phat[2] * (-logSF) ** (1. / phat[0]) - elif ix == 2: - phati = (x - phat[1]) / (-logSF) ** (1. / phat[0]) - else: - raise IndexError('Index to the fixed parameter is out of bounds') - return phati - def _pdf(self, x, c): return c*pow(x, c-1)*exp(-pow(x, c)) @@ -1762,33 +1742,6 @@ class genpareto_gen(rv_continuous): %(example)s """ - - def _link(self, x, logSF, phat, ix): - # Reference - # Stuart Coles (2004) - # "An introduction to statistical modelling of extreme values". - # Springer series in statistics - c, loc, scale = phat - if ix == 2: - # Reorganizing w.r.t.scale, Eq. 4.13 and 4.14, pp 81 in - # Coles (2004) gives - # link = -(x-loc)*c/expm1(-c*logSF) - if c != 0.0: - phati = (x - loc) * c / expm1(-c * logSF) - else: - phati = -(x - loc) / logSF - elif ix == 1: - if c != 0: - phati = x + scale * expm1(c * logSF) / c - else: - phati = x + scale * logSF - elif ix == 0: - raise NotImplementedError( - 'link(x,logSF,phat,i) where i=0 is not implemented!') - else: - raise IndexError('Index to the fixed parameter is out of bounds') - return phati - def _argcheck(self, c): c = asarray(c) self.b = _lazywhere(c < 0, (c,), @@ -1832,57 +1785,6 @@ class genpareto_gen(rv_continuous): scale = m * ((m / s) ** 2 + 1) / 2 return shape, loc, scale - def hessian_nnlf(self, theta, x, eps=None): - 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) - cond0 = (x <= self.a) | (x >= self.b) - if any(cond0): - np = self.numargs + 2 - return valarray((np, np), value=nan) - eps = _EPS - c = args[0] - n = len(x) - if abs(c) > eps: - cx = c * x - sumlog1pcx = sum(log1p(cx)) - # LL = n*log(scale) + (1-1/k)*sumlog1mkxn - r = x / (1.0 + cx) - sumix = sum(1.0 / (1.0 + cx) ** 2.0) - - sumr = sum(r) - sumr2 = sum(r ** 2.0) - H11 = -2 * sumlog1pcx / c ** 3 + 2 * \ - sumr / c ** 2 + (1.0 + 1.0 / c) * sumr2 - H22 = c * (c + 1) * sumix / scale ** 2.0 - H33 = (n - 2 * (c + 1) * sumr + - c * (c + 1) * sumr2) / scale ** 2.0 - H12 = -sum((1 - x) / ((1 + cx) ** 2.0)) / scale - H23 = -(c + 1) * sumix / scale ** 2.0 - H13 = -(sumr - (c + 1) * sumr2) / scale - - else: # c == 0 - sumx = sum(x) - # LL = n*log(scale) + sumx; - - sumx2 = sum(x ** 2.0) - H11 = -(2 / 3) * sum(x ** 3.0) + sumx2 - H22 = 0.0 - H12 = -(n - sum(x)) / scale - H23 = -n * 1.0 / scale ** 2.0 - H33 = (n - 2 * sumx) / scale ** 2.0 - H13 = -(sumx - sumx2) / scale - - # Hessian matrix - H = [[H11, H12, H13], [H12, H22, H23], [H13, H23, H33]] - return asarray(H) - def __stats(self, c): # return None,None,None,None k = -c @@ -1950,25 +1852,9 @@ class genexpon_gen(rv_continuous): %(example)s """ - - def _link(self, x, logSF, phat, ix): - _a, b, c, loc, scale = phat - xn = (x - loc) / scale - fact1 = (xn + expm1(-c * xn) / c) - if ix == 0: - phati = b * fact1 + logSF - elif ix == 1: - phati = (phat[0] - logSF) / fact1 - elif ix in [2, 3, 4]: - raise NotImplementedError('Only implemented for index in [0,1]!') - else: - raise IndexError('Index to the fixed parameter is out of bounds') - - return phati - def _pdf(self, x, a, b, c): - return (a + b*(-special.expm1(-c*x)))*exp((-a-b)*x + - b*(-special.expm1(-c*x))/c) + return (a + b*(-special.expm1(-c*x))) * exp((-a-b)*x + + b*(-special.expm1(-c*x))/c) def _cdf(self, x, a, b, c): return -special.expm1((-a-b)*x + b*(-special.expm1(-c*x))/c) @@ -4248,15 +4134,6 @@ class rayleigh_gen(rv_continuous): %(example)s """ - - def _link(self, x, logSF, phat, ix): - if ix == 1: - return x - phat[0] / sqrt(-2.0 * logSF) - elif ix == 0: - return x - phat[1] * sqrt(-2.0 * logSF) - else: - raise IndexError('Index to the fixed parameter is out of bounds') - def _rvs(self): return chi.rvs(2, size=self._size, random_state=self._random_state) @@ -4314,18 +4191,6 @@ class truncrayleigh_gen(rv_continuous): def _argcheck(self, c): return (c >= 0) - def _link(self, x, logSF, phat, ix): - c, loc, scale = phat - if ix == 2: - return x - loc / (sqrt(c * c - 2 * logSF) - c) - elif ix == 1: - return x - scale * (sqrt(c * c - 2 * logSF) - c) - elif ix == 0: - xn = (x - loc) / scale - return - 2 * logSF / xn - xn / 2.0 - else: - raise IndexError('Index to the fixed parameter is out of bounds') - def _fitstart(self, data, args=None): if args is None: args = (0.0,) * self.numargs @@ -5167,7 +5032,7 @@ _distn_names, _distn_gen_names = get_distribution_names(pairs, rv_continuous) __all__ = _distn_names + _distn_gen_names -if __name__=='__main__': +if __name__ == '__main__': v = genextreme.logpdf(np.inf, 0) v2 = genextreme.logpdf(-np.inf, 0) v2 = genextreme.logpdf(-100, 0) diff --git a/wafo/stats/estimation.py b/wafo/stats/estimation.py index 6fbd769..fb220ef 100644 --- a/wafo/stats/estimation.py +++ b/wafo/stats/estimation.py @@ -21,19 +21,115 @@ from scipy.linalg import pinv2 from scipy import optimize import numpy as np -from numpy import (alltrue, arange, zeros, log, sqrt, exp, +from numpy import (alltrue, arange, zeros, log, sqrt, exp, expm1, atleast_1d, any, asarray, nan, pi, isfinite) from numpy import flatnonzero as nonzero __all__ = ['Profile', 'FitDistribution'] + floatinfo = np.finfo(float) arr = asarray all = alltrue # @ReservedAssignment +def _expon_link(x, logSF, phat, ix): + if ix == 1: + return - (x - phat[0]) / logSF + elif ix == 0: + return x + phat[1] * logSF + else: + raise IndexError('Index to the fixed parameter is out of bounds') + + +def _weibull_min_link(x, logSF, phat, ix): + if ix == 0: + phati = log(-logSF) / log((x - phat[1]) / phat[2]) + elif ix == 1: + phati = x - phat[2] * (-logSF) ** (1. / phat[0]) + elif ix == 2: + phati = (x - phat[1]) / (-logSF) ** (1. / phat[0]) + else: + raise IndexError('Index to the fixed parameter is out of bounds') + return phati + + +def _genpareto_link(x, logSF, phat, ix): + # Reference + # Stuart Coles (2004) + # "An introduction to statistical modelling of extreme values". + # Springer series in statistics + c, loc, scale = phat + if ix == 2: + # Reorganizing w.r.t.scale, Eq. 4.13 and 4.14, pp 81 in + # Coles (2004) gives + # link = -(x-loc)*c/expm1(-c*logSF) + if c != 0.0: + phati = (x - loc) * c / expm1(-c * logSF) + else: + phati = -(x - loc) / logSF + elif ix == 1: + if c != 0: + phati = x + scale * expm1(c * logSF) / c + else: + phati = x + scale * logSF + elif ix == 0: + raise NotImplementedError( + 'link(x,logSF,phat,i) where i=0 is not implemented!') + else: + raise IndexError('Index to the fixed parameter is out of bounds') + return phati + + +def _genexpon_link(x, logSF, phat, ix): + _a, b, c, loc, scale = phat + xn = (x - loc) / scale + fact1 = (xn + expm1(-c * xn) / c) + if ix == 0: + phati = b * fact1 + logSF + elif ix == 1: + phati = (phat[0] - logSF) / fact1 + elif ix in [2, 3, 4]: + raise NotImplementedError('Only implemented for index in [0,1]!') + else: + raise IndexError('Index to the fixed parameter is out of bounds') + + return phati + + +def _rayleigh_link(x, logSF, phat, ix): + if ix == 1: + return x - phat[0] / sqrt(-2.0 * logSF) + elif ix == 0: + return x - phat[1] * sqrt(-2.0 * logSF) + else: + raise IndexError('Index to the fixed parameter is out of bounds') + + +def _trunclayleigh_link(x, logSF, phat, ix): + c, loc, scale = phat + if ix == 2: + return x - loc / (sqrt(c * c - 2 * logSF) - c) + elif ix == 1: + return x - scale * (sqrt(c * c - 2 * logSF) - c) + elif ix == 0: + xn = (x - loc) / scale + return - 2 * logSF / xn - xn / 2.0 + else: + raise IndexError('Index to the fixed parameter is out of bounds') + + +LINKS = dict(expon=_expon_link, + weibull_min=_weibull_min_link, + frechet_r=_weibull_min_link, + genpareto=_genpareto_link, + genexpon=_genexpon_link, + rayleigh=_rayleigh_link, + trunclayleigh=_trunclayleigh_link) + + def chi2isf(p, df): return special.chdtri(df, p) @@ -195,7 +291,7 @@ class Profile(object): self.profile_par = not (self.profile_x or self.profile_logSF) if self.link is None: - self.link = self.fit_dist.dist.link + self.link = LINKS[self.fit_dist.dist.name] if self.profile_par: self._local_link = self._par_link self.xlabel = 'phat(%d)' % self.i_fixed @@ -604,13 +700,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 +893,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]