[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)
         #putmask(cx,cx<-1,-1.0)
         logpdf = where((cx==inf) | (cx==-1),-inf,-(x+cx)*log1pxdx(cx))
         putmask(logpdf,(c==-1) & (x==1.0),0.0)
         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.


More information about the Scipy-tickets mailing list