From db8e835b51aee0d641c6cc1e0a85b05f8acf2d39 Mon Sep 17 00:00:00 2001 From: "Per.Andreas.Brodtkorb" Date: Wed, 6 Oct 2010 12:11:22 +0000 Subject: [PATCH] Fixed bugs in: link method of frechet_r_gen class _reduce_func method of rv_continuous and FitDistribution classes _myprbfun method in Profile class --- pywafo/src/wafo/stats/distributions.py | 15 ++++---- pywafo/src/wafo/stats/estimation.py | 51 ++++++++++++++------------ 2 files changed, 35 insertions(+), 31 deletions(-) diff --git a/pywafo/src/wafo/stats/distributions.py b/pywafo/src/wafo/stats/distributions.py index 047c937..02040ef 100644 --- a/pywafo/src/wafo/stats/distributions.py +++ b/pywafo/src/wafo/stats/distributions.py @@ -2098,10 +2098,10 @@ class rv_continuous(rv_generic): def _reduce_func(self, args, kwds): args = list(args) - Nargs = len(args) - 2 + Nargs = len(args) fixedn = [] - index = range(Nargs) + [-2, -1] - names = ['f%d' % n for n in range(Nargs)] + ['floc', 'fscale'] + index = range(Nargs) + names = ['f%d' % n for n in range(Nargs-2)] + ['floc', 'fscale'] x0 = args[:] for n, key in zip(index, names): if kwds.has_key(key): @@ -3192,16 +3192,15 @@ class frechet_r_gen(rv_continuous): def link(self, x, logSF, phat, ix): u = phat[1] if ix == 0: - phati = (x - phat[1]) / (-logSF) ** (1. / phat[2]) - elif ix == 2: - phati = log(-logSF) / log((x - phat[1]) / phat[0]) + phati = log(-logSF) / log((x - phat[1]) / phat[2]) elif ix == 1: - phati = x - phat[0] * (-logSF) ** (1. / phat[2]); + 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)) def _logpdf(self, x, c): diff --git a/pywafo/src/wafo/stats/estimation.py b/pywafo/src/wafo/stats/estimation.py index 86893f2..3d59ccf 100644 --- a/pywafo/src/wafo/stats/estimation.py +++ b/pywafo/src/wafo/stats/estimation.py @@ -133,7 +133,7 @@ class rv_frozen(object): class Profile(object): ''' Profile Log- likelihood or Product Spacing-function. which can be used for constructing confidence interval for - either phat(i), probability or quantile. + either phat[i], probability or quantile. Parameters ---------- @@ -290,17 +290,18 @@ class Profile(object): self.xlabel = self.xlabel + ' (' + fit_dist.dist.name + ')' - ## Check that par are actually at the optimum + phatfree = phatv[self.i_free].copy() - foundNewphat = False + mylogfun = self._nlogfun if True: - phatfree = optimize.fmin(mylogfun, phatfree, args=(phatv[self.i_fixed],) , disp=0) - LLt = -mylogfun(phatfree,phatv[self.i_fixed]) + ## Check that par are actually at the optimum + phatfree = optimize.fmin(mylogfun, phatfree, args=(p_opt,) , disp=0) + LLt = -mylogfun(phatfree, p_opt) if LLt>Lmax: foundNewphat = True - warnings.warn('Something wrong with fit') + warnings.warn('The fitted parameters does not provide the optimum fit. Something wrong with fit') dL = Lmax-LLt self.alpha_cross_level -= dL self.Lmax = LLt @@ -363,7 +364,6 @@ class Profile(object): pcov = self.fit_dist.par_cov[i_notfixed, :][:, i_notfixed] pvar = sum(numpy.dot(drl, pcov) * drl) - #p_crit = _WAFODIST.norm.isf(self.alpha / 2.0) * sqrt(numpy.ravel(pvar)) * 1.5 p_crit = -norm_ppf(self.alpha / 2.0) * sqrt(numpy.ravel(pvar)) * 1.5 if self.pmin == None: self.pmin = p_opt - 5.0 * p_crit @@ -390,7 +390,7 @@ class Profile(object): def _myprbfun(self, phatnotfixed): mphat = self._par.copy() mphat[self.i_notfixed] = phatnotfixed; - return self.fit_dist.dist.sf(self.x, *mphat) + return self.fit_dist.dist.logsf(self.x, *mphat) def _nlogfun(self, free_par, fix_par): @@ -535,7 +535,7 @@ class FitDistribution(rv_frozen): >>> x = phat.isf(SF) # CI for x - >>> Lx = Profile(phat, i=0,x=x,link=phat.dist.link) + >>> Lx = phat.profile(i=0,x=x,link=phat.dist.link) >>> Lx.plot() >>> x_ci = Lx.get_bounds(alpha=0.2) @@ -641,10 +641,10 @@ class FitDistribution(rv_frozen): def _reduce_func(self, args, kwds): args = list(args) - Nargs = len(args) - 2 + Nargs = len(args) fixedn = [] - index = range(Nargs) + [-2, -1] - names = ['f%d' % n for n in range(Nargs)] + ['floc', 'fscale'] + index = range(Nargs) + names = ['f%d' % n for n in range(Nargs-2)] + ['floc', 'fscale'] x0 = args[:] for n, key in zip(index, names): if kwds.has_key(key): @@ -706,16 +706,8 @@ class FitDistribution(rv_frozen): optimizer = getattr(optimize, optimizer) except AttributeError: raise ValueError, "%s is not a valid optimizer" % optimizer - loglike = numpy.inf - for num_times in range(5): - vals = optimizer(func,x0,args=(ravel(data),),disp=0) - ll = func(vals, ravel(data)) - if ll