WGEVFIT Parameter estimates for GEV data CALL: [phat,cov] = wgevfit(data,method,start,plotflag) phat = Estimated parameters [k s m] = [shape scale location] (see wgevcdf) cov = asymptotic covariance matrix of estimates data = an one-dimensional data set method = a string, describing the method of estimation 'PWM' = Probability Weighted Moments (default) 'ML' = Maximum Likelihood estimates start = starting values for 'ML' method (default 'PWM' estimate) plotflag = 0, do not plot > 0, plot the empiricial distribution function and the estimated cdf (see empdistr for options)(default) The variances of the ML estimates are usually smaller than those of the PWM estimates. However, it is recommended to first use the PWM method since it works for a wider range of parameter values. If method = 'ml' and start is missing, then pwm is used to give the starting values start. Example: R = wgevrnd(0.2,2,7.5,200,1); [phat,cov] = wgevfit(R,'pwm',[],0) [phat,cov] = wgevfit(R,'ml',[],0) See also wgevcdf, empdistr, wgevlike
Computes and plots the empirical CDF | |
Generalized Extreme Value cumulative distribution function | |
Parameter estimates for GEV data | |
Inverse of the Normal distribution function | |
Remove trailing blanks. | |
Display message and abort function. | |
Gamma function. | |
Matrix inverse. | |
Convert string to lowercase. | |
Average or mean value. | |
Create/alter OPTIM OPTIONS structure. | |
Compare strings ignoring case. | |
Graph title. |
% CHAPTER5 contains the commands used in Chapter 5 of the tutorial | |
Parameter estimates for GEV data | |
Tests whether the shape parameter in a GEV is equal to zero |
001 function [phat,cov,pci] = wgevfitwgevfit(data,method,start,plotflag) 002 %WGEVFIT Parameter estimates for GEV data 003 % 004 % CALL: [phat,cov] = wgevfit(data,method,start,plotflag) 005 % 006 % phat = Estimated parameters 007 % [k s m] = [shape scale location] (see wgevcdf) 008 % cov = asymptotic covariance matrix of estimates 009 % data = an one-dimensional data set 010 % method = a string, describing the method of estimation 011 % 'PWM' = Probability Weighted Moments (default) 012 % 'ML' = Maximum Likelihood estimates 013 % start = starting values for 'ML' method 014 % (default 'PWM' estimate) 015 % plotflag = 0, do not plot 016 % > 0, plot the empiricial distribution function and the 017 % estimated cdf (see empdistr for options)(default) 018 % 019 % The variances of the ML estimates are usually smaller than those of 020 % the PWM estimates. However, it is recommended to first use the PWM 021 % method since it works for a wider range of parameter values. 022 % If method = 'ml' and start is missing, then pwm is used to 023 % give the starting values start. 024 % 025 % Example: 026 % R = wgevrnd(0.2,2,7.5,200,1); 027 % [phat,cov] = wgevfit(R,'pwm',[],0) 028 % [phat,cov] = wgevfit(R,'ml',[],0) 029 % 030 % See also wgevcdf, empdistr, wgevlike 031 032 % References 033 % 034 % Prescott, P. and Walden, A.T. (1980) 035 % Maximum likelihood estimation of the parameters of the generalized 036 % extreme-value distribution 037 % Biometrika (67), pp. 723-724 038 % 039 % Hosking, J.R.M, Wallis, J.R. and Wood E.F. (1985) 040 % Estimation of the generalized extreme-value distribution by the 041 % method of probability-weighted moments 042 % Technometrics (27), pp. 251-261 043 % 044 % Borg, S. (1992) 045 % XS - a statistical program package in Splus for extreme-value 046 % analysis. 047 % Dept. of Mathematical Statistics, Lund University, 1992:E2 048 049 % Tested on: Matlab 5.3 050 % History: 051 % Revised by pab 13.06.2001 052 % - Replaced 053 % [phat,dummy,Converged] = fminsearch('wgevlike',mlstart,....); 054 % with 055 % [phat,dummy,Converged] = feval('fminsearch','wgevlike',mlstart,...); 056 % to avoid "Unknown function referenced: fminsearch" error on matlab 5.2 057 % and below. 058 % 059 % Revised by PJ 02-Apr-2001 060 % Method 'ml' now works with new format of fminsearch. 061 % Revised by jr 30.09.1999 062 % Modified by PJ 08-Mar-2000 063 % Added 'hold off' in plotting. 064 % Added routine 'gevll'. Now method 'ml' works. 065 % revised ms 14.06.2000 066 % - updated header info 067 % - changed name to wgevfit (from gevfit) 068 % - added w* to used WAFO-files 069 % - enabled consistent use of 3 and 4 arguments 070 % - enabled use of empty start in ML 071 % revised pab 29.10.2000 072 % - added nargchk 073 074 error(nargchk(1,4,nargin)) 075 if (nargin < 2)|isempty(method), method = 'pwm'; end 076 if (nargin < 4)|isempty(plotflag), plotflag = 1; end 077 078 data = data(:)'; % make sure it is a vector 079 cov=[]; 080 pci=[]; 081 switch lower(method), 082 case 'pwm', 083 w=[ -0.4, 1.6637, 1.3355, 1.1405, 1.8461, 1.1628, 2.9092; 084 -0.3, 1.4153, 0.8912, 0.5640, 1.2574, 0.4442, 1.4090; 085 -0.2, 1.3322, 0.6727, 0.3926, 1.0013, 0.2697, 0.9139; 086 -0.1, 1.2915, 0.5104, 0.3245, 0.8440, 0.2240, 0.6815; 087 0.0, 1.2686, 0.3704, 0.2992, 0.7390, 0.2247, 0.5633; 088 0.1, 1.2551, 0.2411, 0.2966, 0.6708, 0.2447, 0.5103; 089 0.2, 1.2474, 0.1177, 0.3081, 0.6330, 0.2728, 0.5021; 090 0.3, 1.2438, -0.0023, 0.3297, 0.6223, 0.3033, 0.5294; 091 0.4, 1.2433, -0.1205, 0.3592, 0.6368, 0.3329, 0.5880]; 092 093 n = length(data); 094 sortdata = sort(data); 095 koeff1 = ((1:n)-1)/(n-1); 096 koeff2 = koeff1.*((1:n)-2)/(n-2); 097 b2 = (koeff2*sortdata')/n; 098 b1 = (koeff1*sortdata')/n; 099 b0 = mean(sortdata); 100 z = (2*b1-b0)/(3*b2-b0)-log(2)/log(3); 101 shape = 7.8590*z+2.9554*z^2; 102 scale = (2*b1-b0)*shape/(gamma(1+shape)*(1-2^(-shape))); 103 location = b0+scale*(gamma(1+shape)-1)/shape; 104 105 phat=[shape scale location]; 106 % if (abs(shape)>=0.5) changed by ms to deal with shapes in (.4,.5) 107 if (abs(shape)>=0.4) 108 disp(' The estimate of shape is not within the range where ') 109 disp(' the PWM estimator is valid.') 110 % elseif (abs(shape)<=0.5) changed by ms to deal with shapes in (.4,.5) 111 elseif (abs(shape)<=0.4) & nargout>1 112 % calculate covariance matrix of estimates with 113 % linear interpolation between rows 114 i1 = sum((w(:,1)<=shape)); 115 i2 = 10-sum((w(:,1)>=shape)); 116 w_s = w(i1,:)+(shape-w(i1,1))*(w(i2,:)-w(i1,:))/(w(i2,1)-w(i1,1)); 117 W=[w_s(7),w_s(6),w_s(4);w_s(6),w_s(5),w_s(3);w_s(4),w_s(3),w_s(2)]; 118 Cov=[1,scale,scale;scale,scale^2,scale^2;scale,scale^2,scale^2]; 119 cov=1/n*W.*Cov; 120 end 121 122 case 'ml', 123 % ml.gev <- function(data, shape=NA, scale=NA, location=NA) { 124 % J. Statist. Comput. Simul. 16, 241-250 125 if (nargin<3)|isempty(start) % compute starting values by pwm 126 mlstart=wgevfit(data,'pwm',[],0);% Added ms 127 else, 128 mlstart=start; 129 end 130 131 options(1) = 0; % Display off 132 options(2) = 1.e-5; % Termination tolerance on the function value 1.e-5]; 133 options(3) = 1.e-5; % Termination tolerance on X 134 options(14)= 500; % Maximum number of iterations allowed 135 136 % Solve the ML-equation 137 138 try 139 % For Matlab 5.3 and higher ??? 140 [options1] = optimset('disp','off','TolX',options(2),'TolFun',options(3),'MaxIter',options(14)); 141 [phat,dummy,Converged] = feval('fminsearch','wgevlike',mlstart,options1,data); 142 143 catch 144 % For Matlab 5.2 and lower ??? 145 [phat,option] = feval('fmins','wgevlike',mlstart,options,[],data); 146 Converged = (option(14) <= options(14)); 147 end 148 149 150 if ~Converged %(options(10)==500), 151 disp(' ML-routine did not converge in 500 steps.'); cov=[]; 152 elseif nargout>1 153 % Calculate the covariance matrix by inverting the observed information. 154 shape = phat(1); 155 scale = phat(2); 156 location = phat(3); 157 y = -log(1-shape*(data-location)/scale)/shape; 158 P = length(data)-sum(exp(-y)); 159 Q = sum(exp((shape-1)*y)+(shape-1)*exp(shape*y)); 160 R = length(data)-sum(y.*(1-exp(-y))); 161 dLda = -(P+Q)/scale/shape; 162 dLdb = -Q/scale; 163 dLdk = -(R-(P+Q)/shape)/shape; 164 dPdb = -sum(exp((shape-1)*y))/scale; 165 dPda = sum(exp(-y))/scale/shape+dPdb/shape; 166 dQdb = sum(exp((2*shape-1)*y)); 167 dQdb = (dQdb+shape*sum(exp(2*shape*y)))*(1-shape)/scale; 168 dQda = sum(exp((shape-1)*y)); 169 dQda = -(dQda+shape*sum(shape*y))*(1-shape)/scale/shape + dQdb/shape; 170 dRda = length(data)-sum(exp(shape*y))+sum(y.*exp(-y))-sum(exp(-y)); 171 dRda = dRda+sum(exp((shape-1)*y))-sum(y.*exp((shape-1)*y)); 172 dRda = -dRda/scale/shape; 173 dRdk = (sum(y)-sum(y.*exp(-y))+sum(y.*y.*exp(-y))-scale*dRda)/shape; 174 dRdb = (sum(exp(shape*y).*(1-exp(-y)+y.*exp(-y))))/scale; 175 d2Lda2 = -dLda/scale-(dPda+dQda)/scale/shape; 176 d2Ldab = -(dPdb+dQdb)/scale/shape; 177 d2Ldak = -(dRda+dLda+scale*d2Lda2)/shape; 178 d2Ldb2 = -dQdb/scale; 179 d2Ldbk = -(dRdb+scale*d2Ldab)/shape; 180 d2Ldk2 = -(dRdk+dLdk+scale*d2Ldak)/shape; 181 V = - [d2Ldk2, d2Ldak, d2Ldbk; 182 d2Ldak, d2Lda2, d2Ldab; 183 d2Ldbk, d2Ldab, d2Ldb2]; 184 cov = real(inv(V)); 185 end 186 end 187 188 189 if (nargout>2)&~isempty(cov) 190 alpha2 = ones(1,3)*0.05/2; 191 var = diag(cov).'; 192 pci = wnorminv([alpha2;1-alpha2], [phat;phat],[var;var]); 193 end 194 195 % Plotting, JR 9909 196 if plotflag 197 sd=sort(data); 198 empdistr(sd(:),[sd(:), wgevcdf(sd(:),shape,scale,location)],plotflag) 199 if strcmpi(method,'pwm'), 200 str = 'PWM'; 201 else 202 str = 'ML'; 203 end 204 title([deblank(['Empirical and GEV estimated cdf (',str,' method)'])]) 205 end 206 207 208 209 210
Comments or corrections to the WAFO group