From 629ed411c966b883573b80fddc0273ed69419b85 Mon Sep 17 00:00:00 2001 From: "Per.Andreas.Brodtkorb" Date: Thu, 13 Mar 2014 12:46:32 +0000 Subject: [PATCH] pep8 + some refactoring of XXXX.link() methods... --- pywafo/src/wafo/stats/_continuous_distns.py | 291 ++++-------------- .../src/wafo/stats/_distn_infrastructure.py | 47 ++- pywafo/src/wafo/stats/estimation.py | 46 +-- 3 files changed, 124 insertions(+), 260 deletions(-) diff --git a/pywafo/src/wafo/stats/_continuous_distns.py b/pywafo/src/wafo/stats/_continuous_distns.py index 427d4c2..79c822f 100644 --- a/pywafo/src/wafo/stats/_continuous_distns.py +++ b/pywafo/src/wafo/stats/_continuous_distns.py @@ -36,10 +36,10 @@ from ._tukeylambda_stats import (tukeylambda_variance as _tlvar, tukeylambda_kurtosis as _tlkurt) from ._distn_infrastructure import ( - rv_continuous, valarray, - _skew, _kurtosis, _lazywhere, - _ncx2_log_pdf, _ncx2_pdf, _ncx2_cdf, -) + rv_continuous, valarray, + _skew, _kurtosis, _lazywhere, + _ncx2_log_pdf, _ncx2_pdf, _ncx2_cdf, + ) from ._constants import _XMIN, _EULER, _ZETA3, _EPS @@ -68,13 +68,11 @@ __all__ = [ # Kolmogorov-Smirnov one-sided and two-sided test statistics class ksone_gen(rv_continuous): - """General Kolmogorov-Smirnov one-sided test. %(default)s """ - def _cdf(self, x, n): return 1.0 - special.smirnov(n, x) @@ -84,13 +82,11 @@ ksone = ksone_gen(a=0.0, name='ksone') class kstwobign_gen(rv_continuous): - """Kolmogorov-Smirnov two-sided test for large N. %(default)s """ - def _cdf(self, x): return 1.0 - special.kolmogorov(x) @@ -144,7 +140,6 @@ def _norm_isf(q): class norm_gen(rv_continuous): - """A normal continuous random variable. The location (loc) keyword specifies the mean. @@ -161,7 +156,6 @@ class norm_gen(rv_continuous): %(example)s """ - def _rvs(self): return mtrand.standard_normal(self._size) @@ -230,7 +224,6 @@ norm = norm_gen(name='norm') class alpha_gen(rv_continuous): - """An alpha continuous random variable. %(before_notes)s @@ -246,7 +239,6 @@ class alpha_gen(rv_continuous): %(example)s """ - def _pdf(self, x, a): return 1.0 / (x ** 2) / special.ndtr(a) * _norm_pdf(a - 1.0 / x) @@ -265,7 +257,6 @@ alpha = alpha_gen(a=0.0, name='alpha') class anglit_gen(rv_continuous): - """An anglit continuous random variable. %(before_notes)s @@ -281,7 +272,6 @@ class anglit_gen(rv_continuous): %(example)s """ - def _pdf(self, x): return cos(2 * x) @@ -301,7 +291,6 @@ anglit = anglit_gen(a=-pi / 4, b=pi / 4, name='anglit') class arcsine_gen(rv_continuous): - """An arcsine continuous random variable. %(before_notes)s @@ -316,7 +305,6 @@ class arcsine_gen(rv_continuous): %(example)s """ - def _pdf(self, x): return 1.0 / pi / sqrt(x * (1 - x)) @@ -342,7 +330,6 @@ class FitDataError(ValueError): # This exception is raised by, for example, beta_gen.fit when both floc # and fscale are fixed and there are values in the data not in the open # interval (floc, floc+fscale). - def __init__(self, distr, lower, upper): self.args = ( "Invalid values in `data`. Maximum likelihood " @@ -355,7 +342,6 @@ class FitDataError(ValueError): class FitSolverError(RuntimeError): # This exception is raised by, for example, beta_gen.fit when # optimize.fsolve returns with ier != 1. - def __init__(self, mesg): emsg = "Solver for the MLE equations failed to converge: " emsg += mesg.replace('\n', '') @@ -386,7 +372,6 @@ def _beta_mle_ab(theta, n, s1, s2): class beta_gen(rv_continuous): - """A beta continuous random variable. %(before_notes)s @@ -403,7 +388,6 @@ class beta_gen(rv_continuous): %(example)s """ - def _rvs(self, a, b): return mtrand.beta(a, b, self._size) @@ -546,7 +530,6 @@ beta = beta_gen(a=0.0, b=1.0, name='beta') class betaprime_gen(rv_continuous): - """A beta prime continuous random variable. %(before_notes)s @@ -563,7 +546,6 @@ class betaprime_gen(rv_continuous): %(example)s """ - def _rvs(self, a, b): u1 = gamma.rvs(a, size=self._size) u2 = gamma.rvs(b, size=self._size) @@ -602,7 +584,6 @@ betaprime = betaprime_gen(a=0.0, name='betaprime') class bradford_gen(rv_continuous): - """A Bradford continuous random variable. %(before_notes)s @@ -618,7 +599,6 @@ class bradford_gen(rv_continuous): %(example)s """ - def _pdf(self, x, c): return c / (c * x + 1.0) / log1p(c) @@ -662,7 +642,6 @@ bradford = bradford_gen(a=0.0, b=1.0, name='bradford') class burr_gen(rv_continuous): - """A Burr continuous random variable. %(before_notes)s @@ -682,7 +661,6 @@ class burr_gen(rv_continuous): %(example)s """ - def _pdf(self, x, c, d): return c * d * (x ** (-c - 1.0)) * ((1 + x ** (-c * 1.)) ** (-d - 1.)) @@ -699,7 +677,6 @@ burr = burr_gen(a=0.0, name='burr') class fisk_gen(burr_gen): - """A Fisk continuous random variable. The Fisk distribution is also known as the log-logistic distribution, and @@ -714,7 +691,6 @@ class fisk_gen(burr_gen): %(example)s """ - def _pdf(self, x, c): return burr_gen._pdf(self, x, c, 1.0) @@ -734,7 +710,6 @@ fisk = fisk_gen(a=0.0, name='fisk') # median = loc class cauchy_gen(rv_continuous): - """A Cauchy continuous random variable. %(before_notes)s @@ -748,7 +723,6 @@ class cauchy_gen(rv_continuous): %(example)s """ - def _pdf(self, x): return 1.0 / pi / (1.0 + x * x) @@ -776,7 +750,6 @@ cauchy = cauchy_gen(name='cauchy') class chi_gen(rv_continuous): - """A chi continuous random variable. %(before_notes)s @@ -831,7 +804,6 @@ chi = chi_gen(a=0.0, name='chi') # Chi-squared (gamma-distributed with loc=0 and scale=2 and shape=df/2) class chi2_gen(rv_continuous): - """A chi-squared continuous random variable. %(before_notes)s @@ -845,7 +817,6 @@ class chi2_gen(rv_continuous): %(example)s """ - def _rvs(self, df): return mtrand.chisquare(df, self._size) @@ -885,7 +856,6 @@ chi2 = chi2_gen(a=0.0, name='chi2') class cosine_gen(rv_continuous): - """A cosine continuous random variable. %(before_notes)s @@ -902,7 +872,6 @@ class cosine_gen(rv_continuous): %(example)s """ - def _pdf(self, x): return 1.0 / 2 / pi * (1 + cos(x)) @@ -919,7 +888,6 @@ cosine = cosine_gen(a=-pi, b=pi, name='cosine') class dgamma_gen(rv_continuous): - """A double gamma continuous random variable. %(before_notes)s @@ -935,7 +903,6 @@ class dgamma_gen(rv_continuous): %(example)s """ - def _rvs(self, a): u = mtrand.random_sample(size=self._size) return (gamma.rvs(a, size=self._size) * where(u >= 0.5, 1, -1)) @@ -967,7 +934,6 @@ dgamma = dgamma_gen(name='dgamma') class dweibull_gen(rv_continuous): - """A double Weibull continuous random variable. %(before_notes)s @@ -981,7 +947,6 @@ class dweibull_gen(rv_continuous): %(example)s """ - def _rvs(self, c): u = mtrand.random_sample(size=self._size) return weibull_min.rvs(c, size=self._size) * (where(u >= 0.5, 1, -1)) @@ -1017,7 +982,6 @@ dweibull = dweibull_gen(name='dweibull') # Exponential (gamma distributed with a=1.0, loc=loc and scale=scale) class expon_gen(rv_continuous): - """An exponential continuous random variable. %(before_notes)s @@ -1038,28 +1002,13 @@ class expon_gen(rv_continuous): """ - def link(self, x, logSF, phat, ix): - ''' Link for x,SF and parameters of Exponential distribution - - CALL phati = expon.link(x,logSF,phat,i) - - phati = parameter i as function of x, logSF and phat(j) where j ~= i - x = quantile - logSF = logarithm of the survival probability - - LINK is a function connecting the quantile (x) and the survival - probability (R) with the fixed distribution parameter, i.e.: - phat(i) = link(x,logSF,phat,i), - where logSF = log(Prob(X>x;phat)). - - Example % See proflog - - See also profile - ''' + 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 mtrand.standard_exponential(self._size) @@ -1094,7 +1043,6 @@ expon = expon_gen(a=0.0, name='expon') class exponweib_gen(rv_continuous): - """An exponentiated Weibull continuous random variable. %(before_notes)s @@ -1111,10 +1059,8 @@ class exponweib_gen(rv_continuous): %(example)s """ - def _pdf(self, x, a, c): - exc = exp(-x ** c) - return a * c * (1 - exc) ** asarray(a - 1) * exc * x ** (c - 1) + return exp(self._logpdf(x, a, c)) def _logpdf(self, x, a, c): exc = exp(-x ** c) @@ -1131,7 +1077,6 @@ exponweib = exponweib_gen(a=0.0, name='exponweib') class exponpow_gen(rv_continuous): - """An exponential power continuous random variable. %(before_notes)s @@ -1153,7 +1098,6 @@ class exponpow_gen(rv_continuous): %(example)s """ - def _pdf(self, x, b): xbm1 = x ** (b - 1.0) xb = xbm1 * x @@ -1161,7 +1105,7 @@ class exponpow_gen(rv_continuous): def _logpdf(self, x, b): xb = x ** (b - 1.0) * x - return 1 + log(b) + (b - 1.0) * log(x) + xb - exp(xb) + return 1 + log(b) + special.xlogy(b - 1.0, x) + xb - exp(xb) def _cdf(self, x, b): return -expm1(-expm1(x ** b)) @@ -1178,7 +1122,6 @@ exponpow = exponpow_gen(a=0.0, name='exponpow') class fatiguelife_gen(rv_continuous): - """A fatigue-life (Birnbaum-Sanders) continuous random variable. %(before_notes)s @@ -1195,7 +1138,6 @@ class fatiguelife_gen(rv_continuous): %(example)s """ - def _rvs(self, c): z = mtrand.standard_normal(self._size) x = 0.5 * c * z @@ -1234,7 +1176,6 @@ fatiguelife = fatiguelife_gen(a=0.0, name='fatiguelife') class foldcauchy_gen(rv_continuous): - """A folded Cauchy continuous random variable. %(before_notes)s @@ -1250,7 +1191,6 @@ class foldcauchy_gen(rv_continuous): %(example)s """ - def _rvs(self, c): return abs(cauchy.rvs(loc=c, size=self._size)) @@ -1266,7 +1206,6 @@ foldcauchy = foldcauchy_gen(a=0.0, name='foldcauchy') class f_gen(rv_continuous): - """An F continuous random variable. %(before_notes)s @@ -1284,7 +1223,6 @@ class f_gen(rv_continuous): %(example)s """ - def _rvs(self, dfn, dfd): return mtrand.f(dfn, dfd, self._size) @@ -1358,7 +1296,6 @@ f = f_gen(a=0.0, name='f') # Half-normal is folded normal with shape-parameter c=0. class foldnorm_gen(rv_continuous): - """A folded normal continuous random variable. %(before_notes)s @@ -1374,7 +1311,6 @@ class foldnorm_gen(rv_continuous): %(example)s """ - def _argcheck(self, c): return (c >= 0) @@ -1413,7 +1349,6 @@ foldnorm = foldnorm_gen(a=0.0, name='foldnorm') # a limiting value distribution. # class frechet_r_gen(rv_continuous): - """A Frechet right (or Weibull minimum) continuous random variable. %(before_notes)s @@ -1435,8 +1370,7 @@ class frechet_r_gen(rv_continuous): """ - def link(self, x, logSF, phat, ix): - #u = phat[1] + def _link(self, x, logSF, phat, ix): if ix == 0: phati = log(-logSF) / log((x - phat[1]) / phat[2]) elif ix == 1: @@ -1451,7 +1385,7 @@ class frechet_r_gen(rv_continuous): return c * pow(x, c - 1) * exp(-pow(x, c)) def _logpdf(self, x, c): - return log(c) + (c - 1) * log(x) - pow(x, c) + return log(c) + special.xlogy(c - 1, x) - pow(x, c) def _cdf(self, x, c): return -expm1(-pow(x, c)) @@ -1476,7 +1410,6 @@ weibull_min = frechet_r_gen(a=0.0, name='weibull_min') class frechet_l_gen(rv_continuous): - """A Frechet left (or Weibull maximum) continuous random variable. %(before_notes)s @@ -1497,7 +1430,6 @@ class frechet_l_gen(rv_continuous): %(example)s """ - def _pdf(self, x, c): return c * pow(-x, c - 1) * exp(-pow(-x, c)) @@ -1528,7 +1460,6 @@ weibull_max = frechet_l_gen(b=0.0, name='weibull_max') class genlogistic_gen(rv_continuous): - """A generalized logistic continuous random variable. %(before_notes)s @@ -1544,7 +1475,6 @@ class genlogistic_gen(rv_continuous): %(example)s """ - def _pdf(self, x, c): return exp(self._logpdf(x, c)) @@ -1580,7 +1510,6 @@ def log1pxdx(x): class genpareto_gen(rv_continuous): - """A generalized Pareto continuous random variable. %(before_notes)s @@ -1602,29 +1531,28 @@ class genpareto_gen(rv_continuous): """ - def link(self, x, logSF, phat, ix): + def _link(self, x, logSF, phat, ix): # Reference # Stuart Coles (2004) # "An introduction to statistical modelling of extreme values". # Springer series in statistics - - u = phat[1] - if ix == 0: - raise ValueError( - 'link(x,logSF,phat,i) where i=0 is not implemented!') - elif ix == 2: - # Reorganizing w.r.t. phat[2] (scale), Eq. 4.13 and 4.14, pp 81 in + 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-phat[1]).*phat[0]/expm1(phat[0]*logSF) - if phat[0] != 0.0: - phati = (x - u) * phat[0] / expm1(-phat[0] * logSF) + # link = -(x-loc)*c/expm1(-c*logSF) + if c != 0.0: + phati = (x - loc) * c / expm1(-c * logSF) else: - phati = -(x - u) / logSF + phati = -(x - loc) / logSF elif ix == 1: - if phat[0] != 0: - phati = x + phat[2] * expm1(phat[0] * logSF) / phat[0] + if c != 0: + phati = x + scale * expm1(c * logSF) / c else: - phati = x + phat(2) * logSF + 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 @@ -1789,7 +1717,6 @@ genpareto = genpareto_gen(a=0.0, name='genpareto') class genexpon_gen(rv_continuous): - """A generalized exponential continuous random variable. %(before_notes)s @@ -1815,17 +1742,19 @@ class genexpon_gen(rv_continuous): """ - def link(self, x, logSF, phat, ix): - xn = (x - phat[3]) / phat[4] - b = phat[1] - c = phat[2] + 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('Only implemented for ix in [1,2]!') + raise IndexError('Index to the fixed parameter is out of bounds') + return phati def _pdf(self, x, a, b, c): @@ -1842,7 +1771,6 @@ genexpon = genexpon_gen(a=0.0, name='genexpon') class genextreme_gen(rv_continuous): - """A generalized extreme value continuous random variable. %(before_notes)s @@ -1863,7 +1791,6 @@ class genextreme_gen(rv_continuous): %(example)s """ - def _argcheck(self, c): min = np.minimum # @ReservedAssignment max = np.maximum # @ReservedAssignment @@ -1903,9 +1830,9 @@ class genextreme_gen(rv_continuous): def _ppf(self, q, c): x = -log(-log(q)) - return where((c == 0) * (x == x), x, -expm1(-c * x) / c) - # return _lazywhere((x==x) & (c != 0), (x, c), - # lambda x, c: -expm1(-c*x) / c, x) + #return where((c == 0) * (x == x), x, -expm1(-c * x) / c) + return _lazywhere((x==x) & (c != 0), (x, c), + lambda x, c: -expm1(-c*x) / c, x) def _stats(self, c): g = lambda n: gam(n * c + 1) @@ -2006,7 +1933,6 @@ def _digammainv(y): # gamma(df/2, 0, 2) is the chi2 distribution with df degrees of freedom. class gamma_gen(rv_continuous): - """A gamma continuous random variable. %(before_notes)s @@ -2040,7 +1966,6 @@ class gamma_gen(rv_continuous): %(example)s """ - def _rvs(self, a): return mtrand.standard_gamma(a, self._size) @@ -2141,7 +2066,6 @@ gamma = gamma_gen(a=0.0, name='gamma') class erlang_gen(gamma_gen): - """An Erlang continuous random variable. %(before_notes)s @@ -2201,7 +2125,6 @@ erlang = erlang_gen(a=0.0, name='erlang') class gengamma_gen(rv_continuous): - """A generalized gamma continuous random variable. %(before_notes)s @@ -2217,12 +2140,14 @@ class gengamma_gen(rv_continuous): %(example)s """ - def _argcheck(self, a, c): return (a > 0) & (c != 0) def _pdf(self, x, a, c): - return abs(c) * exp((c * a - 1) * log(x) - x ** c - gamln(a)) + return exp(self._logpdf(x, a, c)) + + def _logpdf(self, x, a, c): + return log(abs(c)) + special.xlogy(c * a - 1, x) - x ** c - gamln(a) def _cdf(self, x, a, c): val = special.gammainc(a, x ** c) @@ -2289,7 +2214,6 @@ genhalflogistic = genhalflogistic_gen(a=0.0, name='genhalflogistic') class gompertz_gen(rv_continuous): - """A Gompertz (or truncated Gumbel) continuous random variable. %(before_notes)s @@ -2305,7 +2229,6 @@ class gompertz_gen(rv_continuous): %(example)s """ - def _pdf(self, x, c): return exp(self._logpdf(x, c)) @@ -2324,7 +2247,6 @@ gompertz = gompertz_gen(a=0.0, name='gompertz') class gumbel_r_gen(rv_continuous): - """A right-skewed Gumbel continuous random variable. %(before_notes)s @@ -2346,7 +2268,6 @@ class gumbel_r_gen(rv_continuous): %(example)s """ - def _pdf(self, x): return exp(self._logpdf(x)) @@ -2372,7 +2293,6 @@ gumbel_r = gumbel_r_gen(name='gumbel_r') class gumbel_l_gen(rv_continuous): - """A left-skewed Gumbel continuous random variable. %(before_notes)s @@ -2394,7 +2314,6 @@ class gumbel_l_gen(rv_continuous): %(example)s """ - def _pdf(self, x): return exp(self._logpdf(x)) @@ -2417,7 +2336,6 @@ gumbel_l = gumbel_l_gen(name='gumbel_l') class halfcauchy_gen(rv_continuous): - """A Half-Cauchy continuous random variable. %(before_notes)s @@ -2433,7 +2351,6 @@ class halfcauchy_gen(rv_continuous): %(example)s """ - def _pdf(self, x): return 2.0 / pi / (1.0 + x * x) @@ -2455,7 +2372,6 @@ halfcauchy = halfcauchy_gen(a=0.0, name='halfcauchy') class halflogistic_gen(rv_continuous): - """A half-logistic continuous random variable. %(before_notes)s @@ -2471,7 +2387,6 @@ class halflogistic_gen(rv_continuous): %(example)s """ - def _pdf(self, x): return exp(self._logpdf(x)) @@ -2502,7 +2417,6 @@ halflogistic = halflogistic_gen(a=0.0, name='halflogistic') class halfnorm_gen(rv_continuous): - """A half-normal continuous random variable. %(before_notes)s @@ -2520,7 +2434,6 @@ class halfnorm_gen(rv_continuous): %(example)s """ - def _rvs(self): return abs(mtrand.standard_normal(size=self._size)) @@ -2547,7 +2460,6 @@ halfnorm = halfnorm_gen(a=0.0, name='halfnorm') class hypsecant_gen(rv_continuous): - """A hyperbolic secant continuous random variable. %(before_notes)s @@ -2561,7 +2473,6 @@ class hypsecant_gen(rv_continuous): %(example)s """ - def _pdf(self, x): return 1.0 / (pi * cosh(x)) @@ -2580,7 +2491,6 @@ hypsecant = hypsecant_gen(name='hypsecant') class gausshyper_gen(rv_continuous): - """A Gauss hypergeometric continuous random variable. %(before_notes)s @@ -2598,7 +2508,6 @@ class gausshyper_gen(rv_continuous): %(example)s """ - def _argcheck(self, a, b, c, z): return (a > 0) & (b > 0) & (c == c) & (z == z) @@ -2616,7 +2525,6 @@ gausshyper = gausshyper_gen(a=0.0, b=1.0, name='gausshyper') class invgamma_gen(rv_continuous): - """An inverted gamma continuous random variable. %(before_notes)s @@ -2634,7 +2542,6 @@ class invgamma_gen(rv_continuous): %(example)s """ - def _pdf(self, x, a): return exp(self._logpdf(x, a)) @@ -2670,7 +2577,6 @@ invgamma = invgamma_gen(a=0.0, name='invgamma') # scale is gamma from DATAPLOT and B from Regress class invgauss_gen(rv_continuous): - """An inverse Gaussian continuous random variable. %(before_notes)s @@ -2690,7 +2596,6 @@ class invgauss_gen(rv_continuous): %(example)s """ - def _rvs(self, mu): return mtrand.wald(mu, 1.0, size=self._size) @@ -2715,7 +2620,6 @@ invgauss = invgauss_gen(a=0.0, name='invgauss') class invweibull_gen(rv_continuous): - """An inverted Weibull continuous random variable. %(before_notes)s @@ -2736,7 +2640,6 @@ class invweibull_gen(rv_continuous): %(example)s """ - def _pdf(self, x, c): xc1 = np.power(x, -c - 1.0) xc2 = np.power(x, -c) @@ -2759,7 +2662,6 @@ invweibull = invweibull_gen(a=0, name='invweibull') class johnsonsb_gen(rv_continuous): - """A Johnson SB continuous random variable. %(before_notes)s @@ -2779,7 +2681,6 @@ class johnsonsb_gen(rv_continuous): %(example)s """ - def _argcheck(self, a, b): return (b > 0) & (a == a) @@ -2796,7 +2697,6 @@ johnsonsb = johnsonsb_gen(a=0.0, b=1.0, name='johnsonb') class johnsonsu_gen(rv_continuous): - """A Johnson SU continuous random variable. %(before_notes)s @@ -2817,7 +2717,6 @@ class johnsonsu_gen(rv_continuous): %(example)s """ - def _argcheck(self, a, b): return (b > 0) & (a == a) @@ -2835,7 +2734,6 @@ johnsonsu = johnsonsu_gen(name='johnsonsu') class laplace_gen(rv_continuous): - """A Laplace continuous random variable. %(before_notes)s @@ -2849,7 +2747,6 @@ class laplace_gen(rv_continuous): %(example)s """ - def _rvs(self): return mtrand.laplace(0, 1, size=self._size) @@ -2860,7 +2757,7 @@ class laplace_gen(rv_continuous): return where(x > 0, 1.0 - 0.5 * exp(-x), 0.5 * exp(x)) def _ppf(self, q): - return where(q > 0.5, -log(2 * (1 - q)), log(2 * q)) + return where(q > 0.5, -log(2) - log1p(-q), log(2 * q)) def _stats(self): return 0, 2, 0, 3 @@ -2871,7 +2768,6 @@ laplace = laplace_gen(name='laplace') class levy_gen(rv_continuous): - """A Levy continuous random variable. %(before_notes)s @@ -2895,7 +2791,7 @@ class levy_gen(rv_continuous): """ def _pdf(self, x): - return 1 / sqrt(2 * pi * x) / x * exp(-1 / (2 * x)) + return 1 / sqrt(2*pi*x) / x * exp(-1/(2*x)) def _cdf(self, x): return 2 * (1 - _norm_cdf(1 / sqrt(x))) @@ -2910,7 +2806,6 @@ levy = levy_gen(a=0.0, name="levy") class levy_l_gen(rv_continuous): - """A left-skewed Levy continuous random variable. %(before_notes)s @@ -2932,7 +2827,6 @@ class levy_l_gen(rv_continuous): %(example)s """ - def _pdf(self, x): ax = abs(x) return 1 / sqrt(2 * pi * ax) / ax * exp(-1 / (2 * ax)) @@ -2951,7 +2845,6 @@ levy_l = levy_l_gen(b=0.0, name="levy_l") class levy_stable_gen(rv_continuous): - """A Levy-stable continuous random variable. %(before_notes)s @@ -2968,7 +2861,6 @@ class levy_stable_gen(rv_continuous): %(example)s """ - def _rvs(self, alpha, beta): sz = self._size TH = uniform.rvs(loc=-pi / 2.0, scale=pi, size=sz) @@ -3004,7 +2896,6 @@ levy_stable = levy_stable_gen(name='levy_stable') class logistic_gen(rv_continuous): - """A logistic (or Sech-squared) continuous random variable. %(before_notes)s @@ -3020,7 +2911,6 @@ class logistic_gen(rv_continuous): %(example)s """ - def _rvs(self): return mtrand.logistic(size=self._size) @@ -3046,7 +2936,6 @@ logistic = logistic_gen(name='logistic') class loggamma_gen(rv_continuous): - """A log gamma continuous random variable. %(before_notes)s @@ -3062,7 +2951,6 @@ class loggamma_gen(rv_continuous): %(example)s """ - def _rvs(self, c): return log(mtrand.gamma(c, size=self._size)) @@ -3088,7 +2976,6 @@ loggamma = loggamma_gen(name='loggamma') class loglaplace_gen(rv_continuous): - """A log-Laplace continuous random variable. %(before_notes)s @@ -3110,7 +2997,6 @@ class loglaplace_gen(rv_continuous): %(example)s """ - def _pdf(self, x, c): cd2 = c / 2.0 c = where(x < 1, c, -c) @@ -3137,7 +3023,6 @@ def _lognorm_logpdf(x, s): class lognorm_gen(rv_continuous): - """A lognormal continuous random variable. %(before_notes)s @@ -3157,7 +3042,6 @@ class lognorm_gen(rv_continuous): %(example)s """ - def _rvs(self, s): return exp(s * mtrand.standard_normal(self._size)) @@ -3195,7 +3079,6 @@ lognorm = lognorm_gen(a=0.0, name='lognorm') class gilbrat_gen(rv_continuous): - """A Gilbrat continuous random variable. %(before_notes)s @@ -3211,7 +3094,6 @@ class gilbrat_gen(rv_continuous): %(example)s """ - def _rvs(self): return exp(mtrand.standard_normal(self._size)) @@ -3246,7 +3128,6 @@ gilbrat = gilbrat_gen(a=0.0, name='gilbrat') class maxwell_gen(rv_continuous): - """A Maxwell continuous random variable. %(before_notes)s @@ -3269,7 +3150,6 @@ class maxwell_gen(rv_continuous): %(example)s """ - def _rvs(self): return chi.rvs(3.0, size=self._size) @@ -3295,7 +3175,6 @@ maxwell = maxwell_gen(a=0.0, name='maxwell') class mielke_gen(rv_continuous): - """A Mielke's Beta-Kappa continuous random variable. %(before_notes)s @@ -3311,7 +3190,6 @@ class mielke_gen(rv_continuous): %(example)s """ - def _pdf(self, x, k, s): return k * x ** (k - 1.0) / (1.0 + x ** s) ** (1.0 + k * 1.0 / s) @@ -3325,7 +3203,6 @@ mielke = mielke_gen(a=0.0, name='mielke') class nakagami_gen(rv_continuous): - """A Nakagami continuous random variable. %(before_notes)s @@ -3342,7 +3219,6 @@ class nakagami_gen(rv_continuous): %(example)s """ - def _pdf(self, x, nu): return (2 * nu ** nu / gam(nu) * (x ** (2 * nu - 1.0)) * exp(-nu * x * x)) @@ -3364,7 +3240,6 @@ nakagami = nakagami_gen(a=0.0, name="nakagami") class ncx2_gen(rv_continuous): - """A non-central chi-squared continuous random variable. %(before_notes)s @@ -3381,7 +3256,6 @@ class ncx2_gen(rv_continuous): %(example)s """ - def _rvs(self, df, nc): return mtrand.noncentral_chisquare(df, nc, self._size) @@ -3413,7 +3287,6 @@ ncx2 = ncx2_gen(a=0.0, name='ncx2') class ncf_gen(rv_continuous): - """A non-central F distribution continuous random variable. %(before_notes)s @@ -3434,7 +3307,6 @@ class ncf_gen(rv_continuous): %(example)s """ - def _rvs(self, dfn, dfd, nc): return mtrand.noncentral_f(dfn, dfd, nc, self._size) @@ -3475,7 +3347,6 @@ ncf = ncf_gen(a=0.0, name='ncf') class t_gen(rv_continuous): - """A Student's T continuous random variable. %(before_notes)s @@ -3493,7 +3364,6 @@ class t_gen(rv_continuous): %(example)s """ - def _rvs(self, df): return mtrand.standard_t(df, size=self._size) @@ -3530,7 +3400,6 @@ t = t_gen(name='t') class nct_gen(rv_continuous): - """A non-central Student's T continuous random variable. %(before_notes)s @@ -3548,7 +3417,6 @@ class nct_gen(rv_continuous): %(example)s """ - def _argcheck(self, df, nc): return (df > 0) & (nc == nc) @@ -3616,7 +3484,6 @@ nct = nct_gen(name="nct") class pareto_gen(rv_continuous): - """A Pareto continuous random variable. %(before_notes)s @@ -3632,7 +3499,6 @@ class pareto_gen(rv_continuous): %(example)s """ - def _pdf(self, x, b): return b * x ** (-b - 1) @@ -3675,7 +3541,6 @@ pareto = pareto_gen(a=1.0, name="pareto") class lomax_gen(rv_continuous): - """A Lomax (Pareto of the second kind) continuous random variable. %(before_notes)s @@ -3694,7 +3559,6 @@ class lomax_gen(rv_continuous): %(example)s """ - def _pdf(self, x, c): return c * 1.0 / (1.0 + x) ** (c + 1.0) @@ -3723,7 +3587,6 @@ lomax = lomax_gen(a=0.0, name="lomax") class pearson3_gen(rv_continuous): - """A pearson type III continuous random variable. %(before_notes)s @@ -3756,7 +3619,6 @@ class pearson3_gen(rv_continuous): Aviation Loads Data", Office of Aviation Research (2003). """ - def _preprocess(self, x, skew): # The real 'loc' and 'scale' are handled in the calling pdf(...). The # local variables 'loc' and 'scale' within pearson3._pdf are set to @@ -3778,7 +3640,7 @@ class pearson3_gen(rv_continuous): invmask = ~mask beta = 2.0 / (skew[invmask] * scale) - alpha = (scale * beta) ** 2 + alpha = (scale * beta)**2 zeta = loc - alpha / beta transx = beta * (x[invmask] - zeta) @@ -3792,11 +3654,11 @@ class pearson3_gen(rv_continuous): return np.ones(np.shape(skew), dtype=bool) def _stats(self, skew): - _, _x, _transx, skew, _mask, _invmask, beta, alpha, zeta = ( + ans, x, transx, skew, mask, invmask, beta, alpha, zeta = ( self._preprocess([1], skew)) m = zeta + alpha / beta - v = alpha / (beta ** 2) - s = 2.0 / (alpha ** 0.5) * np.sign(beta) + v = alpha / (beta**2) + s = 2.0 / (alpha**0.5) * np.sign(beta) k = 6.0 / alpha return m, v, s, k @@ -3825,7 +3687,7 @@ class pearson3_gen(rv_continuous): return ans def _cdf(self, x, skew): - ans, x, transx, skew, mask, invmask, _beta, alpha, _zeta = ( + ans, x, transx, skew, mask, invmask, beta, alpha, _zeta = ( self._preprocess(x, skew)) ans[mask] = _norm_cdf(x[mask]) @@ -3837,22 +3699,21 @@ class pearson3_gen(rv_continuous): self._preprocess([0], skew)) if mask[0]: return mtrand.standard_normal(self._size) - ans = mtrand.standard_gamma(alpha, self._size) / beta + zeta + ans = mtrand.standard_gamma(alpha, self._size)/beta + zeta if ans.size == 1: return ans[0] return ans def _ppf(self, q, skew): - ans, q, _transq, skew, mask, invmask, beta, alpha, zeta = ( + ans, q, transq, skew, mask, invmask, beta, alpha, zeta = ( self._preprocess(q, skew)) ans[mask] = _norm_ppf(q[mask]) - ans[invmask] = special.gammaincinv(alpha, q[invmask]) / beta + zeta + ans[invmask] = special.gammaincinv(alpha, q[invmask])/beta + zeta return ans pearson3 = pearson3_gen(name="pearson3") class powerlaw_gen(rv_continuous): - """A power-function continuous random variable. %(before_notes)s @@ -3870,12 +3731,11 @@ class powerlaw_gen(rv_continuous): %(example)s """ - def _pdf(self, x, a): return a * x ** (a - 1.0) def _logpdf(self, x, a): - return log(a) + (a - 1) * log(x) + return log(a) + special.xlogy(a - 1, x) def _cdf(self, x, a): return x ** (a * 1.0) @@ -3898,7 +3758,6 @@ powerlaw = powerlaw_gen(a=0.0, b=1.0, name="powerlaw") class powerlognorm_gen(rv_continuous): - """A power log-normal continuous random variable. %(before_notes)s @@ -3916,7 +3775,6 @@ class powerlognorm_gen(rv_continuous): %(example)s """ - def _pdf(self, x, c, s): return (c / (x * s) * _norm_pdf(log(x) / s) * pow(_norm_cdf(-log(x) / s), c * 1.0 - 1.0)) @@ -3930,7 +3788,6 @@ powerlognorm = powerlognorm_gen(a=0.0, name="powerlognorm") class powernorm_gen(rv_continuous): - """A power normal continuous random variable. %(before_notes)s @@ -3947,7 +3804,6 @@ class powernorm_gen(rv_continuous): %(example)s """ - def _pdf(self, x, c): return (c * _norm_pdf(x) * (_norm_cdf(-x) ** (c - 1.0))) @@ -3963,7 +3819,6 @@ powernorm = powernorm_gen(name='powernorm') class rdist_gen(rv_continuous): - """An R-distributed continuous random variable. %(before_notes)s @@ -3979,7 +3834,6 @@ class rdist_gen(rv_continuous): %(example)s """ - def _pdf(self, x, c): return (np.power((1.0 - x ** 2), c / 2.0 - 1) / special.beta(0.5, c / 2.0)) @@ -4001,7 +3855,6 @@ rdist = rdist_gen(a=-1.0, b=1.0, name="rdist") class rayleigh_gen(rv_continuous): - """A Rayleigh continuous random variable. %(before_notes)s @@ -4020,12 +3873,13 @@ class rayleigh_gen(rv_continuous): """ - def link(self, x, logSF, phat, ix): - rv_continuous.link.__doc__ + def _link(self, x, logSF, phat, ix): if ix == 1: return x - phat[0] / sqrt(-2.0 * logSF) - else: + 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) @@ -4077,16 +3931,17 @@ 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] + def _link(self, x, logSF, phat, ix): + c, loc, scale = phat if ix == 2: - return x - phat[1] / (sqrt(c * c - 2 * logSF) - c) + return x - loc / (sqrt(c * c - 2 * logSF) - c) elif ix == 1: - return x - phat[2] * (sqrt(c * c - 2 * logSF) - c) + return x - scale * (sqrt(c * c - 2 * logSF) - c) elif ix == 0: - xn = (x - phat[1]) / phat[2] + 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: @@ -4132,7 +3987,6 @@ truncrayleigh = truncrayleigh_gen(a=0.0, name="truncrayleigh", shapes='c') class reciprocal_gen(rv_continuous): - """A reciprocal continuous random variable. %(before_notes)s @@ -4148,7 +4002,6 @@ class reciprocal_gen(rv_continuous): %(example)s """ - def _argcheck(self, a, b): self.a = a self.b = b @@ -4177,7 +4030,6 @@ reciprocal = reciprocal_gen(name="reciprocal") # FIXME: PPF does not work. class rice_gen(rv_continuous): - """A Rice continuous random variable. %(before_notes)s @@ -4193,7 +4045,6 @@ class rice_gen(rv_continuous): %(example)s """ - def _argcheck(self, b): return b >= 0 @@ -4217,7 +4068,6 @@ rice = rice_gen(a=0.0, name="rice") # FIXME: PPF does not work. class recipinvgauss_gen(rv_continuous): - """A reciprocal inverse Gaussian continuous random variable. %(before_notes)s @@ -4233,7 +4083,6 @@ class recipinvgauss_gen(rv_continuous): %(example)s """ - def _rvs(self, mu): return 1.0 / mtrand.wald(mu, 1.0, size=self._size) @@ -4255,7 +4104,6 @@ recipinvgauss = recipinvgauss_gen(a=0.0, name='recipinvgauss') class semicircular_gen(rv_continuous): - """A semicircular continuous random variable. %(before_notes)s @@ -4271,7 +4119,6 @@ class semicircular_gen(rv_continuous): %(example)s """ - def _pdf(self, x): return 2.0 / pi * sqrt(1 - x * x) @@ -4287,7 +4134,6 @@ semicircular = semicircular_gen(a=-1.0, b=1.0, name="semicircular") class triang_gen(rv_continuous): - """A triangular continuous random variable. %(before_notes)s @@ -4305,7 +4151,6 @@ class triang_gen(rv_continuous): %(example)s """ - def _rvs(self, c): return mtrand.triangular(0, c, 1, self._size) @@ -4334,7 +4179,6 @@ triang = triang_gen(a=0.0, b=1.0, name="triang") class truncexpon_gen(rv_continuous): - """A truncated exponential continuous random variable. %(before_notes)s @@ -4350,7 +4194,6 @@ class truncexpon_gen(rv_continuous): %(example)s """ - def _argcheck(self, b): self.b = b return (b > 0) @@ -4386,7 +4229,6 @@ truncexpon = truncexpon_gen(a=0.0, name='truncexpon') class truncnorm_gen(rv_continuous): - """A truncated normal continuous random variable. %(before_notes)s @@ -4403,7 +4245,6 @@ class truncnorm_gen(rv_continuous): %(example)s """ - def _argcheck(self, a, b): self.a = a self.b = b @@ -4445,7 +4286,6 @@ truncnorm = truncnorm_gen(name='truncnorm') # FIXME: RVS does not work. class tukeylambda_gen(rv_continuous): - """A Tukey-Lamdba continuous random variable. %(before_notes)s @@ -4464,7 +4304,6 @@ class tukeylambda_gen(rv_continuous): %(example)s """ - def _argcheck(self, lam): return np.ones(np.shape(lam), dtype=bool) @@ -4494,7 +4333,6 @@ tukeylambda = tukeylambda_gen(name='tukeylambda') class uniform_gen(rv_continuous): - """A uniform continuous random variable. This distribution is constant between `loc` and ``loc + scale``. @@ -4504,7 +4342,6 @@ class uniform_gen(rv_continuous): %(example)s """ - def _rvs(self): return mtrand.uniform(0.0, 1.0, self._size) @@ -4526,7 +4363,6 @@ uniform = uniform_gen(a=0.0, b=1.0, name='uniform') class vonmises_gen(rv_continuous): - """A Von Mises continuous random variable. %(before_notes)s @@ -4550,7 +4386,6 @@ class vonmises_gen(rv_continuous): %(example)s """ - def _rvs(self, kappa): return mtrand.vonmises(0.0, kappa, size=self._size) @@ -4567,7 +4402,6 @@ vonmises_line = vonmises_gen(a=-np.pi, b=np.pi, name='vonmises_line') class wald_gen(invgauss_gen): - """A Wald continuous random variable. %(before_notes)s @@ -4584,7 +4418,6 @@ class wald_gen(invgauss_gen): %(example)s """ - def _rvs(self): return mtrand.wald(1.0, 1.0, size=self._size) @@ -4603,7 +4436,6 @@ wald = wald_gen(a=0.0, name="wald") class wrapcauchy_gen(rv_continuous): - """A wrapped Cauchy continuous random variable. %(before_notes)s @@ -4619,7 +4451,6 @@ class wrapcauchy_gen(rv_continuous): %(example)s """ - def _argcheck(self, c): return (c > 0) & (c < 1) diff --git a/pywafo/src/wafo/stats/_distn_infrastructure.py b/pywafo/src/wafo/stats/_distn_infrastructure.py index f03709c..b9c6224 100644 --- a/pywafo/src/wafo/stats/_distn_infrastructure.py +++ b/pywafo/src/wafo/stats/_distn_infrastructure.py @@ -23,7 +23,7 @@ from scipy import optimize from scipy import integrate # to approximate the pdf of a continuous distribution given its cdf -from scipy.misc import comb, derivative +from scipy.misc import comb, derivative # @UnresolvedImport from numpy import (arange, putmask, ravel, take, ones, sum, shape, product, reshape, zeros, floor, logical_and, log, sqrt, exp, @@ -1553,7 +1553,7 @@ class rv_continuous(rv_generic): for item in ['callparams', 'default', 'before_notes']: tempdict[item] = tempdict[item].replace( "\n%(shapes)s : array_like\n shape parameters", "") - for i in range(2): + for _i in range(2): if self.shapes is None: # necessary because we use %(shapes)s in two forms (w w/o ", ") self.__doc__ = self.__doc__.replace("%(shapes)s, ", "") @@ -1953,13 +1953,40 @@ class rv_continuous(rv_generic): return output[()] return output - def link(self, x, logSF, theta, i): - ''' Return dist. par. no. i as function of quantile (x) and log survival probability (sf) - where - theta is the list containing all parameters including location and scale. + def link(self, x, logSF, theta, i): + ''' + Return theta[i] as function of quantile, survival probability and + theta[j] for j!=i. + + Parameters + ---------- + x : quantile + logSF : logarithm of the survival probability + theta : list + all distribution parameters including location and scale. + + Returns + ------- + theta[i] : real scalar + fixed distribution parameter theta[i] as function of x, logSF and + theta[j] where j != i. + + LINK is a function connecting the fixed distribution parameter theta[i] + with the quantile (x) and the survival probability (SF) and the + remaining free distribution parameters theta[j] for j!=i, i.e.: + theta[i] = link(x, logSF, theta, i), + where logSF = log(Prob(X>x; theta)). + + See also + estimation.Profile ''' - raise ValueError('Link function not implemented for the %s distribution' % self.name) - return None + return self._link(x, logSF, theta, i) + + def _link(self, x, logSF, theta, i): + msg = ('Link function not implemented for the %s distribution' % + self.name) + raise NotImplementedError(msg) + def _nnlf(self, x, *args): return -sum(self._logpdf(x, *args), axis=0) @@ -3030,7 +3057,7 @@ class rv_discrete(rv_generic): for item in ['callparams', 'default', 'before_notes']: tempdict[item] = tempdict[item].replace( "\n%(shapes)s : array_like\n shape parameters", "") - for i in range(2): + for _i in range(2): if self.shapes is None: # necessary because we use %(shapes)s in two forms (w w/o ", ") self.__doc__ = self.__doc__.replace("%(shapes)s, ", "") @@ -3498,7 +3525,7 @@ class rv_discrete(rv_generic): else: invfac = 1.0 - tot = 0.0 + #tot = 0.0 low, upp = self._ppf(0.001, *args), self._ppf(0.999, *args) low = max(min(-suppnmin, low), lb) upp = min(max(suppnmin, upp), ub) diff --git a/pywafo/src/wafo/stats/estimation.py b/pywafo/src/wafo/stats/estimation.py index 8a29501..808856c 100644 --- a/pywafo/src/wafo/stats/estimation.py +++ b/pywafo/src/wafo/stats/estimation.py @@ -152,8 +152,8 @@ class Profile(object): with ML or MPS estimated distribution parameters. **kwds : named arguments with keys i : scalar integer - defining which distribution parameter to profile, i.e. which - parameter to keep fixed (default first non-fixed parameter) + defining which distribution parameter to keep fixed in the + profiling process (default first non-fixed parameter) pmin, pmax : real scalars Interval for either the parameter, phat(i), prb, or x, used in the optimization of the profile function (default is based on the @@ -166,8 +166,8 @@ class Profile(object): log survival probability,i.e., SF = Prob(X>x;phat) (default None) link : function connecting the x-quantile and the survival probability (SF) with the fixed distribution parameter, i.e.: - self.par[i] = link(x,logSF,self.par,i), where - logSF = log(Prob(X>x;phat)). + self.par[i] = link(x, logSF, self.par, i), where + logSF = log(Prob(X>x;phat)). This means that if: 1) x is not None then x is profiled 2) logSF is not None then logSF is profiled @@ -236,22 +236,21 @@ class Profile(object): self.fit_dist = fit_dist self.data = None self.args = None - self.title = 'Profile log' + self.title = '' self.xlabel = '' - self.ylabel = '' + self.ylabel = 'Profile log' (self.i_fixed, self.N, self.alpha, self.pmin, self.pmax, self.x, self.logSF, self.link) = map( kwds.get, - ['i', 'N', 'alpha', 'pmin', - 'pmax', 'x', 'logSF', 'link'], + ['i', 'N', 'alpha', 'pmin', 'pmax', 'x', 'logSF', 'link'], [i0, 100, 0.05, None, None, None, None, None]) - self.ylabel = '%g%s CI' % (100 * (1.0 - self.alpha), '%') + self.title = '%g%s CI' % (100 * (1.0 - self.alpha), '%') if fit_dist.method.startswith('ml'): - self.title = self.title + 'likelihood' + self.ylabel = self.ylabel + 'likelihood' Lmax = fit_dist.LLmax elif fit_dist.method.startswith('mps'): - self.title = self.title + ' product spacing' + self.ylabel = self.ylabel + ' product spacing' Lmax = fit_dist.LPSmax else: raise ValueError( @@ -267,7 +266,7 @@ class Profile(object): self.i_fixed = atleast_1d(self.i_fixed) if 1 - isnotfixed[self.i_fixed]: - raise ValueError( + raise IndexError( "Index i must be equal to an index to one of the free " + "parameters.") @@ -295,7 +294,7 @@ class Profile(object): self.xlabel = 'phat(%d)' % self.i_fixed p_opt = self._par[self.i_fixed] elif self.profile_x: - self.logSF = log(fit_dist.sf(self.x)) + self.logSF = fit_dist.logsf(self.x) self._local_link = lambda fix_par, par: self.link( fix_par, self.logSF, par, self.i_fixed) self.xlabel = 'x' @@ -522,16 +521,23 @@ class Profile(object): CI = ecross(self.args, self.data, ind[[0, -1]], cross_level) return CI - def plot(self): + def plot(self, axis=None): ''' Plot profile function with 100(1-alpha)% CI ''' - plotbackend.plot( + if axis is None: + axis = plotbackend.gca() + + p_ci = self.get_bounds(self.alpha) + axis.plot( self.args, self.data, - self.args[[0, -1]], [self.Lmax, ] * 2, 'r', - self.args[[0, -1]], [self.alpha_cross_level, ] * 2, 'r') - plotbackend.title(self.title) - plotbackend.ylabel(self.ylabel) - plotbackend.xlabel(self.xlabel) + self.args[[0, -1]], [self.Lmax, ] * 2, 'r--', + self.args[[0, -1]], [self.alpha_cross_level, ] * 2, 'r--') + axis.vlines(p_ci, ymin=axis.get_ylim()[0], + ymax=self.Lmax, #self.alpha_cross_level, + color='r', linestyles='--') + axis.set_title(self.title) + axis.set_ylabel(self.ylabel) + axis.set_xlabel(self.xlabel) def _discretize_adaptive(fun, a, b, tol=0.005, n=5):