WGEVCDF Generalized Extreme Value cumulative distribution function CALL: F = wgevcdf(x,k,s,m); F = distribution function evaluated at x k = shape parameter in the GEV s = scale parameter in the GEV, s>0 (default 1) m = location parameter in the GEV (default 0) The Generalized Extreme Value distribution is defined by its cdf exp( - (1 - k(x-m)/s))^1/k) ), k~=0 F(x;k,s,m) = exp( -exp(-(x-m)/s) ), k==0 for x>s/k+m (when k<=0) and x<m+s/k (when k>0). Example: x = linspace(0,15,200); p1 = wgevcdf(x,0.8,1,11); p2 = wgevcdf(x,0.8,2,11); p3 = wgevcdf(x,0.5,1,11); p4 = wgevcdf(x,0.5,2,11); plot(x,p1,x,p2,x,p3,x,p4)
Check if all input arguments are either scalar or of common size. | |
Display message and abort function. | |
Not-a-Number. |
% CHAPTER5 contains the commands used in Chapter 5 of the tutorial | |
Parameter estimates for GEV data |
001 function F = wgevcdf(x,k,s,m); 002 %WGEVCDF Generalized Extreme Value cumulative distribution function 003 % 004 % CALL: F = wgevcdf(x,k,s,m); 005 % 006 % F = distribution function evaluated at x 007 % k = shape parameter in the GEV 008 % s = scale parameter in the GEV, s>0 (default 1) 009 % m = location parameter in the GEV (default 0) 010 % 011 % The Generalized Extreme Value distribution is defined by its cdf 012 % 013 % exp( - (1 - k(x-m)/s))^1/k) ), k~=0 014 % F(x;k,s,m) = 015 % exp( -exp(-(x-m)/s) ), k==0 016 % 017 % for x>s/k+m (when k<=0) and x<m+s/k (when k>0). 018 % 019 % Example: 020 % x = linspace(0,15,200); 021 % p1 = wgevcdf(x,0.8,1,11); p2 = wgevcdf(x,0.8,2,11); 022 % p3 = wgevcdf(x,0.5,1,11); p4 = wgevcdf(x,0.5,2,11); 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 by jr 22.12.1999 032 % revised ms 14.06.2000 033 % - updated header info 034 % - changed name to wgevcdf (from gevcdf) 035 % revised pab 24.10.2000 036 % - added nargchk, comnsize and default values for m, s 037 % revised jr 14.08.2001 038 % - a bug in the last if-statement condition fixed 039 % (thanks to D Eddelbuettel <edd@debian.org>) 040 041 error(nargchk(2,4,nargin)) 042 043 if nargin<4|isempty(m), m=0;end 044 if nargin<3|isempty(s), s=1;end 045 046 [errorcode x k s,m] = comnsize(x,k,s,m); 047 if errorcode > 0 048 error('x, k, s and m must be of common size or scalar.'); 049 end 050 051 epsilon=1e-4; % treshold defining k to zero 052 053 F = zeros(size(x)); 054 %k0 = find(x>=m & abs(k)<=epsilon & s>0); 055 k0 = find( abs(k)<=epsilon & s>0); 056 if any(k0), 057 F(k0) = exp(-exp(-(x(k0)-m(k0))./s(k0))); 058 end 059 060 %k1=find(x>m&(k.*x<s+k.*m)&(abs(k)>epsilon)); 061 k1=find((k.*x<s+k.*m)&(abs(k)>epsilon)); 062 if any(k1), 063 tmp = (1-k(k1).*(x(k1)-m(k1))./s(k1)); 064 F(k1)=exp(-tmp.^(1./k(k1))); 065 end 066 067 k2=find((k.*x>=s+k.*m)&(k>epsilon)); 068 if any(k2), 069 F(k2)=ones(size(k2)); 070 end 071 072 k3 = find(s<=0 ); 073 if any(k3), 074 tmp = NaN; 075 F(k3) = tmp(ones(size(k3))); 076 end 077 return 078 079 % old call 080 epsilon=1e-4; 081 if abs(k) < epsilon, 082 p = exp(-exp(-(x-m)/s)); 083 else 084 if k>0 085 p=ones(size(x)); 086 p=p.*(x>=s/k+m)+(x<s/k+m).*... 087 exp(-(1-k*(x-m)/s).^(1/k)); 088 else 089 p=zeros(size(x)); 090 p=(x>s/k+m).*... 091 exp(-(1-k*(x-m)/s).^(1/k)); 092 end 093 end 094
Comments or corrections to the WAFO group