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:
 empdistr Computes and plots the empirical CDF wgevcdf Generalized Extreme Value cumulative distribution function wgevfit Parameter estimates for GEV data wnorminv Inverse of the Normal distribution function deblank Remove trailing blanks. error Display message and abort function. gamma Gamma function. inv Matrix inverse. lower Convert string to lowercase. mean Average or mean value. optimset Create/alter OPTIM OPTIONS structure. strcmpi Compare strings ignoring case. title Graph title.
This function is called by:
 Chapter5 % CHAPTER5 contains the commands used in Chapter 5 of the tutorial wgevfit Parameter estimates for GEV data wgumbtest Tests whether the shape parameter in a GEV is equal to zero

## 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 %
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
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
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