Home > wafo > wstats > wgevfit.m

wgevfit

PURPOSE ^

Parameter estimates for GEV data

SYNOPSIS ^

[phat,cov,pci] = wgevfit(data,method,start,plotflag)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

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