From 59186edfd6fc3a71b999a9f87113379331ebaaf0 Mon Sep 17 00:00:00 2001 From: Per A Brodtkorb Date: Wed, 8 Jun 2016 16:00:55 +0200 Subject: [PATCH] Improved accuracy and avoids some runtime warnings. --- wafo/stats/_continuous_distns.py | 39 +++++++++++++++++++++----------- 1 file changed, 26 insertions(+), 13 deletions(-) diff --git a/wafo/stats/_continuous_distns.py b/wafo/stats/_continuous_distns.py index 93ea6de..cdce229 100644 --- a/wafo/stats/_continuous_distns.py +++ b/wafo/stats/_continuous_distns.py @@ -32,12 +32,12 @@ except: from scipy.stats._tukeylambda_stats import (tukeylambda_variance as _tlvar, tukeylambda_kurtosis as _tlkurt) -from ._distn_infrastructure import ( +from wafo.stats._distn_infrastructure import ( rv_continuous, valarray, _skew, _kurtosis, _lazywhere, _ncx2_log_pdf, _ncx2_pdf, _ncx2_cdf, get_distribution_names, ) -from ._constants import _XMIN, _EULER, _ZETA3, _XMAX, _LOGXMAX, _EPS +from wafo.stats._constants import _XMIN, _EULER, _ZETA3, _XMAX, _LOGXMAX, _EPS ## Kolmogorov-Smirnov one-sided and two-sided test statistics @@ -2018,12 +2018,12 @@ class genextreme_gen(rv_continuous): return exp(self._logpdf(x, c)) def _logpdf(self, x, c): - cx = c*x - logex2 = where((c == 0)*(x == x), 0.0, special.log1p(-cx)) + cx = _lazywhere((c != 0)*(x == x), (x, c), lambda x, c: c*x, 0.0) + logex2 = special.log1p(-cx) logpex2 = self._loglogcdf(x, c) - # logpex2 = where((c == 0)*(x == x), -x, logex2/c) pex2 = exp(logpex2) # Handle special cases + putmask(logpex2, (c == 0) & (x == -inf), 0.0) logpdf = where((cx == 1) | (cx == -inf), -inf, -pex2+logpex2-logex2) putmask(logpdf, (c == 1) & (x == 1), 0.0) return logpdf @@ -2032,9 +2032,8 @@ class genextreme_gen(rv_continuous): return exp(self._logcdf(x, c)) def _loglogcdf(self, x, c): - return _lazywhere((c == 0)*(x == x), (x, c), - f=lambda x, c: -x, - f2=lambda x, c: special.log1p(-c*x)/c) + return _lazywhere((c != 0) & (x == x), (x, c), + lambda x, c: special.log1p(-c * x)/c, -x) def _logcdf(self, x, c): return -exp(self._loglogcdf(x, c)) @@ -2049,10 +2048,8 @@ class genextreme_gen(rv_continuous): def _isf(self, q, c): x = -log(-special.log1p(-q)) - result = _lazywhere((c == 0)*(x == x), (x, c), - f=lambda x, c: x, - f2=lambda x, c: -special.expm1(-c*x)/c) - return result + 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) @@ -2560,6 +2557,9 @@ class gumbel_l_gen(rv_continuous): def _logpdf(self, x): return x - exp(x) + def _logsf(self, x): + return -exp(x) + def _sf(self, x): return exp(-exp(x)) @@ -3203,6 +3203,9 @@ class logistic_gen(rv_continuous): def _ppf(self, q): return -log1p(-q) + log(q) + def _isf(self, q): + return log1p(-q) - log(q) + def _stats(self): return 0, pi*pi/3.0, 0, 6.0/5.0 @@ -4262,7 +4265,8 @@ class rayleigh_gen(rv_continuous): def _logpdf(self, r): rr2 = r * r / 2.0 - return where(rr2 == inf, - rr2, log(r) - rr2) + return _lazywhere(rr2 != inf, (r, rr2), lambda r, rr2: log(r) - rr2, + -rr2) def _cdf(self, r): return -special.expm1(-0.5 * r**2) @@ -4273,6 +4277,9 @@ class rayleigh_gen(rv_continuous): def _sf(self, r): return exp(-0.5 * r**2) + def _logsf(self, r): + return -0.5 * r**2 + def _isf(self, q): return sqrt(-2 * log(q)) @@ -5158,3 +5165,9 @@ pairs = list(globals().items()) _distn_names, _distn_gen_names = get_distribution_names(pairs, rv_continuous) __all__ = _distn_names + _distn_gen_names + + +if __name__=='__main__': + v = genextreme.logpdf(np.inf, 0) + v2 = genextreme.logpdf(-np.inf, 0) + v2 = genextreme.logpdf(-100, 0)