Home > wafo > wstats > wgpdcdf.m

wgpdcdf

PURPOSE ^

Generalized Pareto cumulative distribution function

SYNOPSIS ^

F = wgpdcdf(x,k,s,m);

DESCRIPTION ^

 WGPDCDF Generalized Pareto cumulative distribution function 
  
  CALL:  F = wgpdcdf(x,k,s,m); 
   
         F = distribution function evaluated at x 
         k = shape parameter in the GPD 
         s = scale parameter in the GPD    (default 1) 
         m = location parameter in the GPD (default 0) 
   
  The Generalized Pareto distribution is defined by its cdf 
  
                 1 - (1-k(x-m)/s)^1/k,  k~=0 
   F(x;k,s,m) = 
                 1 - exp(-(x-m)/s),  k==0 
  
  for x>=m (when k<=0) and 0 <= x-m < s/k (when k>0), s>0. 
  
  Example:  
    x = linspace(0,2,200); 
    p1 = wgpdcdf(x,1.25,1); p2 = wgpdcdf(x,1,1); 
    p3 = wgpdcdf(x,0.75,1); p4 = wgpdcdf(x,0.5,1); 
    plot(x,p1,x,p2,x,p3,x,p4)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function F = wgpdcdf(x,k,s,m); 
002 %WGPDCDF Generalized Pareto cumulative distribution function 
003 % 
004 % CALL:  F = wgpdcdf(x,k,s,m); 
005 %  
006 %        F = distribution function evaluated at x 
007 %        k = shape parameter in the GPD 
008 %        s = scale parameter in the GPD    (default 1) 
009 %        m = location parameter in the GPD (default 0) 
010 %  
011 % The Generalized Pareto distribution is defined by its cdf 
012 % 
013 %                1 - (1-k(x-m)/s)^1/k,  k~=0 
014 %  F(x;k,s,m) = 
015 %                1 - exp(-(x-m)/s),  k==0 
016 % 
017 % for x>=m (when k<=0) and 0 <= x-m < s/k (when k>0), s>0. 
018 % 
019 % Example:  
020 %   x = linspace(0,2,200); 
021 %   p1 = wgpdcdf(x,1.25,1); p2 = wgpdcdf(x,1,1); 
022 %   p3 = wgpdcdf(x,0.75,1); p4 = wgpdcdf(x,0.5,1); 
023 %   plot(x,p1,x,p2,x,p3,x,p4) 
024   
025 % References 
026 %  Johnson  N.L., Kotz S. and Balakrishnan, N. (1994) 
027 %  Continuous Univariate Distributions, Volume 1. Wiley.  
028  
029 % Tested on: Matlab 5.3 
030 % History:  
031 % Revised pab oct 2005 
032 % - limits m<x<s/k replaced with 0 <= x-m < s/k in help header. 
033 % revised pab June 2005 
034 % - distribution for k==0 was wrong, now fixed 
035 % Revised by jr 22.12.1999 
036 % Modified by PJ 08-Mar-2000 
037 %   Hjälptext 
038 % revised ms 14.06.2000 
039 % - updated header info 
040 % - changed name to wgpdcdf (from gpdcdf) 
041 % revised pab 24.10.2000 
042 % - added  nargchk, comnsize and default values for m, s  
043 % revised jr 14.08.2001 
044 % - a bug in the last if-statement condition fixed 
045 %   (thanks to D Eddelbuettel <edd@debian.org>) 
046  
047 error(nargchk(2,4,nargin)) 
048  
049 if nargin<4|isempty(m), m=0;end 
050 if nargin<3|isempty(s), s=1;end 
051  
052 [errorcode x k s,m] = comnsize(x,k,s,m); 
053 if errorcode > 0 
054   error('x, k, s and m must be of common size or scalar.'); 
055 end 
056    
057 epsilon = 1e-4; % treshold defining k to zero 
058  
059  
060 F  = zeros(size(x)); 
061 xm = x-m; 
062 k0 = find((0 < xm) & (abs(k)<=epsilon) & (0<s)); 
063 if any(k0), 
064   F(k0) = 1-exp(-xm(k0)./s(k0)); 
065 end 
066  
067 k1 = find((0<xm) & (k.*xm < s ) & (abs(k)>epsilon) & (0<s)); 
068 if any(k1), 
069   F(k1) = 1-(1-k(k1).*xm(k1)./s(k1)).^(1./k(k1)); 
070 end 
071  
072 k2 = find((k.*xm>=s) & (k>epsilon)); 
073 if any(k2), 
074   F(k2)=ones(size(k2)); 
075 end 
076  
077 k3 = find(s<=0 ); 
078 if any(k3), 
079    tmp   = NaN; 
080    F(k3) = tmp(ones(size(k3))); 
081  end 
082   
083 return 
084 % old call 
085 if abs(k) < epsilon,   
086     p = (x>0).*(1-exp(-x/s)); 
087   elseif k>0,  
088     p=(x>=s/k)+(x>0).*(x<s/k).*... 
089     (1-(1-k*x/s).^(1/k)); 
090   else  
091     p=zeros(size(x)); 
092     p=(x>0).*(1-(1-k*x/s).^(1/k)); 
093 end; 
094  
095  
096  
097  
098

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003