[Scipy-tickets] [SciPy] #768: Add the limiting exponential distribution to the generalized pareto distribution when shape c=0.

SciPy scipy-tickets@scipy....
Mon Dec 1 09:37:29 CST 2008

```#768: Add the limiting exponential distribution to the generalized pareto
distribution when shape c=0.
-------------------------+--------------------------------------------------
Reporter:  pbrod        |        Owner:  somebody
Type:  enhancement  |       Status:  new
Priority:  normal       |    Milestone:
Component:  scipy.stats  |      Version:
Severity:  normal       |   Resolution:
Keywords:               |
-------------------------+--------------------------------------------------
Comment (by pbrod):

A solution to this issue, which also gives improved accuracy in the upper
tail as well as for small values of the shape parameter, c, is given here:

{{{
class genpareto_gen(rv_continuous):
def _argcheck(self, c):
c = arr(c)
self.b = where(0<=c,inf, 1.0/abs(c))
return where(abs(c)==inf, 0, 1)
def _pdf(self, x, c):
cx = where((c==0) & (x==inf),0.0,c*x).clip(min=-1.0)
logpdf = where((cx==inf) | (cx==-1),-inf,-(x+cx)*log1pxdx(cx))
return exp(logpdf)

#%f = exp(-xn)./s;                   % for  k==0
#%f = (1+k.*xn).^(-1./k-1)/s;        % for  k~=0
#%f = exp((-1./k-1).*log1p(kxn))/s  % for  k~=0
#%f = exp((-xn-kxn).*log1p(kxn)./(kxn))/s  % for any k kxn~=inf
#Px = pow(1+c*x,arr(-1.0-1.0/c))
#return Px
def _chf(self,x,c):
cx = c*x
return where((0.0<x) & (-1.0<=cx) & (c!=0),log1p(cx)/c,x)
def _cdf(self, x, c):
log_sf = -self._chf(x,c)
return -expm1(log_sf)
#return 1.0 - pow(1+c*x,arr(-1.0/c))
def _sf(self,x,c):
log_sf = -self._chf(x,c)
return exp(log_sf)
def _isf2(self,log_sf,c):
return where((c!=0) & (-inf<log_sf),expm1(-c*log_sf)/c,-log_sf)
def _ppf(self, q, c):
log_sf = log1p(-q)
return self._isf2(log_sf,c)
def _isf(self,q,c):
log_sf = log(q)
return self._isf2(log_sf,c)
#vals = 1.0/c * (pow(1-q, -c)-1)
#return vals

def _stats(self,c):
#return None,None,None,None
k = -c
m = where(k<-1.0,inf,1.0/(1+k))
v = where(k<-0.5,nan,1.0/((1+k)**2.0*(1+2*k)))
sk = where(k<-1.0/3,nan,2.*(1-k)*sqrt(1+2.0*k)/(1.0 +3.*k))
#% E(X^r) = s^r*(-k)^-(r+1)*gamma(1+r)*gamma(-1/k-r)/gamma(1-1/k)
#%  = s^r*gamma(1+r)./( (1+k)*(1+2*k).*....*(1+r*k))
#% E[(1-k(X-m0)/s)^r] = 1/(1+k*r)

#%Ex3 = (sk.*sqrt(v)+3*m).*v+m^3
#%Ex3 = 6.*s.^3/((1+k).*(1+2*k).*(1+3*k))
r = 4.0;
Ex4 = gam(1.+r)/((1.+k)*(1.+2.*k)*(1.+3.*k)*(1+4.*k))
m1 = m
ku =
where(k<-1./4,nan,(Ex4-4.*sk*v**(3./2)*m1-6*m1**2.*v-m1**4.)/v**2.-3.0)
return m,v,sk,ku
}}}
where

{{{
def log1pxdx(x):
'''Computes Log(1+x)/x
'''
y = where(x==0,1.0,log1p(x)/x)
return where(x==inf,0.0,y)
}}}

--
Ticket URL: <http://scipy.org/scipy/scipy/ticket/768#comment:1>
SciPy <http://www.scipy.org/>
SciPy is open-source software for mathematics, science, and engineering.
```