WGEVPDF Generalized Extreme Value probability density function CALL: f = wgevpdf(x,k,s,m); f = density 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 = wgevpdf(x,0.8,1,11); p2 = wgevpdf(x,0.8,2,11); p3 = wgevpdf(x,0.5,1,11); p4 = wgevpdf(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. |
001 function f = wgevpdf(x,k,s,m); 002 %WGEVPDF Generalized Extreme Value probability density function 003 % 004 % CALL: f = wgevpdf(x,k,s,m); 005 % 006 % f = density 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 = wgevpdf(x,0.8,1,11); p2 = wgevpdf(x,0.8,2,11); 022 % p3 = wgevpdf(x,0.5,1,11); p4 = wgevpdf(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 jr 14.08.2001 032 % - a bug in the last if-statement condition fixed 033 % (thanks to D Eddelbuettel <edd@debian.org>) 034 % revised pab 24.10.2000 035 % - added nargchk, comnsize and default values for m, s 036 % added ms 14.06.2000 037 038 error(nargchk(2,4,nargin)) 039 040 if nargin<4|isempty(m), m=0;end 041 if nargin<3|isempty(s), s=1;end 042 043 [errorcode x k s,m] = comnsize(x,k,s,m); 044 if errorcode > 0 045 error('x, k, s and m must be of common size or scalar.'); 046 end 047 048 epsilon=1e-4; % treshold defining k to zero 049 050 f = zeros(size(x)); 051 k0 = find(x>=m & abs(k)<=epsilon & s>0); 052 if any(k0), 053 tmp=exp(-(x(k0)-m(k0))./s(k0)); 054 f(k0) = exp(-tmp).*tmp./s(k0); 055 end 056 057 k1=find((k.*x<s+k.*m)&(abs(k)>epsilon)); 058 if any(k1), 059 tmp = (1-k(k1).*(x(k1)-m(k1))./s(k1)); 060 f(k1)=exp(-tmp.^(1./k(k1))).*tmp.^(1./k(k1)-1)./s(k1); 061 end 062 063 %k2=find((k.*x>=s+k.*m)&(k>epsilon)); 064 %if any(k2), 065 % f(k2)=ones(size(k2)); 066 %end 067 068 k3 = find(s<=0 ); 069 if any(k3), 070 tmp = NaN; 071 f(k3) = tmp(ones(size(k3))); 072 end 073 return 074 075 076 077 078
Comments or corrections to the WAFO group