WGEVSTAT Mean and variance for the GEV distribution. CALL: [m,v] = wgevstat(k,s,m0) m, v = the mean and variance, respectively k, s, m0 = parameter of the GEV distribution. (see wgevcdf) Mean (m) and variance (v) for the GEV distribution is m=(m0*k+s)/k-s*gamma(k+1)/k and v=(m0*k+s)^2/(k^2)+(-2*m0*k-2*s)*s*gamma(k+1)/(k^2)... +s^2*gamma(2*k+1)/(k^2)-m^2 Note: mean only exists for k>-1 and variance for k>0 Example: [m,v] = wgevstat(0.5,2,0) [m,v] = wgevstat(-0.5,2,0)
Check if all input arguments are either scalar or of common size. | |
Display message and abort function. | |
Gamma function. | |
Not-a-Number. |
001 function [m,v]= wgevstat(k,s,m0); 002 %WGEVSTAT Mean and variance for the GEV distribution. 003 % 004 % CALL: [m,v] = wgevstat(k,s,m0) 005 % 006 % m, v = the mean and variance, respectively 007 % k, s, m0 = parameter of the GEV distribution. (see wgevcdf) 008 % 009 % Mean (m) and variance (v) for the GEV distribution is 010 % 011 % m=(m0*k+s)/k-s*gamma(k+1)/k and 012 % v=(m0*k+s)^2/(k^2)+(-2*m0*k-2*s)*s*gamma(k+1)/(k^2)... 013 % +s^2*gamma(2*k+1)/(k^2)-m^2 014 % 015 % Note: mean only exists for k>-1 and variance for k>0 016 % 017 % Example: 018 % [m,v] = wgevstat(0.5,2,0) 019 % [m,v] = wgevstat(-0.5,2,0) 020 021 % Tested on; Matlab 5.3 022 % History: 023 % revised pab 24.10.2000 024 % - added comnsize, nargchk + default value for s and m0 025 % added ms 09.08.2000 026 027 error(nargchk(1,3,nargin)) 028 if nargin<2, s=1;end 029 if nargin<3, m0=0;end 030 [errorcode k,s,m0] = comnsize(k,s,m0); 031 if errorcode > 0 032 error('k s and m0 must be of common size or scalar.'); 033 end 034 035 % Initialize m and v to zero. 036 m = zeros(size(k)); 037 v=m; 038 ok = (s > 0 & k>-1 ); 039 k1=find(ok); 040 if any(k1), 041 m(k1) = (m0(k1).*k(k1)+s(k1))./k(k1)-s(k1).*gamma(k(k1)+1)./k(k1); 042 end 043 044 ok1 = (s>0 & k>0 ); 045 k2=find(ok1); 046 if any(k2), 047 v(k2) = (m0(k2).*k(k2)+s(k2)).^2./(k(k2).^2)+... 048 (-2*m0(k2).*k(k2)-2*s(k2)).*s(k2).*gamma(k(k2)+1)./... 049 (k(k2).^2)+s(k2).^2*gamma(2*k(k2)+1)./(k(k2).^2)-m(k2).^2; 050 end 051 052 053 k3 = find(~ok); 054 if any(k3) 055 tmp = NaN; 056 m(k3) = tmp(ones(size(k1))); 057 end 058 k4 = find(~ok1); 059 if any(k4) 060 tmp = NaN; 061 v(k4) = tmp(ones(size(k1))); 062 end 063 064 065 066
Comments or corrections to the WAFO group