0001 function [lcEst,Est,R,MSE] = cmat2extralc(param,F,u,method,plotflag)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100 ni = nargin;
0101 no = nargout;
0102 error(nargchk(3,5,ni));
0103
0104 if ni < 4, method = []; end
0105 if ni < 5, plotflag = []; end
0106
0107
0108 if isempty(method), method = 'gpd,ml'; end
0109 if isempty(plotflag), plotflag = 1; end
0110
0111
0112 n = param(3);
0113 paramD = [1 n n];
0114 lc0 = cmat2lc(paramD,F);
0115 uu = levels(param)';
0116 uuD = levels(paramD)';
0117
0118
0119 u_high = u(2);
0120 i_high = min(uuD(uu>u_high));
0121
0122
0123 u_low = u(1);
0124 i_low = max(uuD(uu<u_low));
0125
0126
0127 F1 = F;
0128 for i = i_high:n
0129 F1(i,:) = 0;
0130 end
0131 for j = 1:i_low
0132 F1(:,j) = 0;
0133 end
0134
0135
0136 lc1 = cmat2lc(paramD,F1);
0137
0138 lc=lc1;
0139
0140 [dummy,I]=max(lc(:,2));
0141 i_LCmax = lc(I(1),1);
0142
0143
0144 if plotflag, subplot(2,1,1), end
0145
0146 lcHigh = lc(i_high:end,:);
0147 [lcEst.D.High,Est.D.High,MSE.High] = extrapolate(lcHigh,method,plotflag,i_high-i_LCmax);
0148
0149 if plotflag, title([method ': Extrapolation for high levels']), end
0150
0151
0152 if plotflag, subplot(2,1,2), end
0153
0154 lcLow = lc(1:i_low,:);
0155 lcLow1 = [n+1-flipud(lcLow(:,1)) flipud(lcLow(:,2))];
0156
0157 [lcEst1,Est1,MSE.Low] = extrapolate(lcLow1,method,plotflag,i_LCmax-i_low);
0158 lcEst.D.Low = [n+1-flipud(lcEst1(:,1)) flipud(lcEst1(:,2))];
0159 Est.D.Low = Est1;
0160
0161 if isfield(Est.D.Low,'UpperLimit');
0162 Est.D.Low.LowerLimit = n+1-Est.D.Low.UpperLimit;
0163 rmfield(Est.D.Low,'UpperLimit');
0164 end
0165
0166 if plotflag, title([method ': Extrapolation for low levels']), end
0167
0168
0169 lcEst.D.lc = lc;
0170 lcEst.D.lc(1:i_low,:) = lcEst.D.Low;
0171 lcEst.D.lc(i_high:end,:) = lcEst.D.High;
0172
0173
0174 lcEst.High = [uu(lcEst.D.High(:,1)) lcEst.D.High(:,2)];
0175 lcEst.Low = [uu(lcEst.D.Low(:,1)) lcEst.D.Low(:,2)];
0176 lcEst.lc = [uu(lcEst.D.lc(:,1)) lcEst.D.lc(:,2)];
0177
0178 if no > 1
0179 Est.method = method;
0180 end
0181
0182
0183 if no > 2
0184 R.lc0 = cmat2lc(param,F);
0185 R.lc1 = cmat2lc(param,F1);
0186 R.F1 = F1;
0187 end
0188
0189
0190
0191
0192 function [lcEst,Est,MSE] = extrapolate(lc,method,plotflag,offset)
0193
0194
0195
0196
0197 iu = lc(1,1);
0198
0199
0200 lc1 = [lc(:,1)-iu lc(:,2)];
0201 lc2 = flipud(lc1);
0202 lc3 = lc2;
0203
0204 method = lower(method);
0205 methodShape = method(1:3);
0206 methodEst = method(5:end);
0207
0208 if strcmp(methodShape,'gpd') & ~strcmp(method,'gpd,ml')
0209
0210 x = [];
0211 for k = 2:length(lc3)
0212 nn = lc3(k,2)-lc3(k-1,2);
0213 if nn > 0
0214 x = [x; lc3(k,1)*ones(nn,1)];
0215 end
0216 end
0217 x = x+.5;
0218 end
0219
0220 if strcmp(methodShape,'exp') | strcmp(method,'gpd,ml') | strcmp(method,'ray,ml')
0221 dN = [];
0222 for k = 2:length(lc3)
0223 nn = lc3(k,2)-lc3(k-1,2);
0224 if nn ~= 0
0225 dN = [dN; lc3(k,1) nn];
0226 end
0227 end
0228 NN = [dN(:,1) cumsum(dN(:,2))];
0229 NN = flipud(NN);
0230 dN = flipud(dN);
0231 end
0232
0233
0234 if strcmp(methodShape,'gpd')
0235
0236 if strcmp(methodEst,'ml')
0237
0238 dNN = [dN(:,1)+0.5 dN(:,2)];
0239 n = NN(1,2);
0240
0241
0242 X1=dNN(1,1); Xn=dNN(end,1); Xmean = sum(dNN(:,1).*dNN(:,2))/sum(dNN(:,2));
0243 Eps = 1e-6/Xmean;
0244 t_L = 2*(X1-Xmean)/X1^2;
0245 t_U = 1/Xn-Eps;
0246 x_L = log(1/Xn - t_L);
0247 x_U = log(1/Xn - t_U);
0248
0249 x=linspace(x_L,x_U,10);
0250 for i=1:length(x)
0251 f(i)=wgpdfit_mld(x(i),dNN);
0252 end
0253
0254 I = find(f(1:end-1).*f(2:end) < 0);
0255
0256
0257 if ~isempty(I)
0258 i = I(1);
0259 x_start = [x(i) x(i+1)];
0260 if exist('optimset') >= 2
0261
0262 x_ML = fzero('wgpdfit_mld',x_start,optimset('disp','off'),dNN);
0263 else
0264
0265 x_ML = fzero('wgpdfit_mld',x_start,[],[],dNN);
0266 end
0267 [f,shape,scale] = wgpdfit_mld(x_ML,dNN);
0268
0269
0270 if (shape >= 1.0),
0271 cov=[];
0272 warning([' The estimate of shape (' num2str(shape) ...
0273 ') is not within the range where the ML estimator is valid (shape<1).'])
0274 elseif (shape >= 0.5),
0275 cov=[];
0276 warning([' The estimate of shape (' num2str(shape) ...
0277 ') is not within the range where the ML estimator is asymptotically normal (shape<1/2).'])
0278 else
0279 Vshape = 1-shape;
0280 Vscale = 2*scale^2;
0281 Covari = scale;
0282 cov = (1-shape)/n*[Vshape Covari; Covari Vscale];
0283 end
0284
0285 else
0286 warning([' ML-estimator does not exist for this data. Using MOM instead.'])
0287
0288 Xvar = sum((dNN(:,1)-Xmean).^2.*dNN(:,2))/(sum(dNN(:,2))-1);
0289 shape = (Xmean^2/Xvar-1)/2;
0290 scale = Xmean*(Xmean^2/Xvar+1)/2;
0291 cov = [];
0292 end
0293 parms = [shape scale];
0294
0295 else
0296 if strcmp(methodEst,'ml0')
0297 [parms,cov] = wgpdfit(x,'ml',0);
0298 else
0299 [parms,cov] = wgpdfit(x,methodEst,0);
0300 end
0301 end
0302
0303 Est.k = parms(1);
0304 Est.s = parms(2);
0305 if Est.k>0,
0306 Est.UpperLimit = iu+Est.s/Est.k;
0307 end
0308 Est.cov = cov;
0309 xF = (0:length(lc1)-1)';
0310 F = wgpdcdf(xF,Est.k,Est.s);
0311
0312 elseif strcmp(methodShape,'exp')
0313
0314 switch methodEst
0315
0316 case 'ml'
0317
0318
0319
0320
0321
0322
0323
0324 n = NN(1,2);
0325 Sx = sum(dN(:,2).*(dN(:,1)+0.5));
0326 parms(1) = 0;
0327 parms(2) = Sx/n;
0328
0329 cov = zeros(2);
0330 cov(2,2) = 1/n*parms(2);
0331
0332 case 'mld'
0333
0334 n = NN(1,2);
0335 Sx = sum(dN(:,2).*(dN(:,1)+0.0));
0336 s1 = 1/log(1+n/Sx);
0337
0338 parms(1) = 0;
0339 parms(2) = s1;
0340
0341 cov = zeros(2);
0342 cov(2,2) = 1/n*parms(2);
0343
0344
0345 case {'ls1','wls1'}
0346
0347 X = NN(:,1);
0348 Y = log(NN(:,2))-log(NN(1,2));
0349 n = length(Y);
0350
0351 if strcmp(methodEst,'ls1')
0352 W = eye(n);
0353 else
0354 W = diag(dN(:,2));
0355 end
0356
0357
0358 th = (X'*W*X)\(X'*W*Y);
0359
0360 parms(1) = 0;
0361 parms(2) = -1/th(1);
0362 cov = [];
0363
0364 case {'ls2','wls2'}
0365
0366 X = [ones(size(NN(:,1))) NN(:,1)];
0367 Y = log(NN(:,2));
0368 n = length(Y);
0369
0370 if strcmp(methodEst,'ls2')
0371 W = eye(n);
0372 else
0373 W = diag(dN(:,2));
0374 end
0375
0376
0377 th = (X'*W*X)\(X'*W*Y);
0378
0379 parms(1) = 0;
0380 parms(2) = -1/th(2);
0381 parms(3) = exp(th(1));
0382 cov = [];
0383
0384 case {'ls3','wls3'}
0385
0386 X = [ones(size(NN(:,1))) NN(:,1) NN(:,1).^2];
0387 Y = log(NN(:,2));
0388 n = length(Y);
0389
0390 if strcmp(methodEst,'ls3')
0391 W = eye(n);
0392 else
0393 W = diag(dN(:,2));
0394 end
0395
0396
0397 th = (X'*W*X)\(X'*W*Y);
0398
0399 parms(1) = 0;
0400 parms(2) = -1/th(2);
0401 parms(3) = exp(th(1));
0402 parms(4) = -1/th(3);
0403 Est.s2 = parms(4);
0404 cov = [];
0405
0406 case {'ls4','wls4'}
0407
0408 X = [NN(:,1) NN(:,1).^2];
0409 Y = log(NN(:,2))-log(NN(1,2));
0410 n = length(Y);
0411
0412 if strcmp(methodEst,'ls4')
0413 W = eye(n);
0414 else
0415 W = diag(dN(:,2));
0416 end
0417
0418
0419 th = (X'*W*X)\(X'*W*Y);
0420
0421 parms(1) = 0;
0422 parms(2) = -1/th(1);
0423 parms(4) = -1/th(2);
0424 Est.s2 = parms(4);
0425 cov = [];
0426
0427 end
0428
0429 Est.k = parms(1);
0430 Est.s = parms(2);
0431 Est.cov = cov;
0432 xF = (0:length(lc1)-1)';
0433 switch methodEst
0434 case {'ml', 'mld', 'ls1', 'wls1', 'ls2', 'wls2'}
0435 F = 1 - exp(-xF/Est.s);
0436 case {'ls3', 'wls3', 'ls4', 'wls4'}
0437 F = 1 - exp(-xF/Est.s-xF.^2/Est.s2);
0438 end
0439
0440 elseif strcmp(methodShape,'ray')
0441
0442 n = NN(1,2);
0443 Sx = sum(dN(:,2).*((dN(:,1)+offset+0.5).^2-(offset)^2));
0444 a=sqrt(Sx/n);
0445
0446 Est.s = a;
0447
0448 xF = (0:length(lc1)-1)';
0449 F = 1 - exp(-((xF+offset).^2-offset^2)/a^2);
0450
0451 else
0452
0453 error(['Unknown method: ' method]);
0454
0455 end
0456
0457
0458 switch methodEst
0459 case {'ls2','wls2','ls3','wls3'}
0460 Est.lcu = parms(3);
0461 otherwise
0462 Est.lcu = lc(1,2);
0463 end
0464
0465 lcEst = [xF+iu Est.lcu*(1-F)];
0466 lcEst(end,2) = 0;
0467
0468
0469 if strcmp(method,'gpd,ml') | strcmp(method,'exp,ml') | strcmp(method,'exp,mld') | strcmp(method,'ray,mld')
0470 MSE = 0;
0471
0472 elseif strcmp(methodShape,'gpd')
0473 n=length(x);
0474 q1 = ((1:n)-1/2)/n;
0475 xq1 = wgpdinv(q1,Est.k,Est.s)';
0476 MSE = sum((x-xq1).^2)/n;
0477 elseif strcmp(method,'ray,ml')
0478 MSE = 0;
0479
0480 else
0481 Lam = sqrt(W);
0482 MSE = (Lam*X*th - Lam*Y)'*(Lam*X*th - Lam*Y)/sum(diag(W));
0483 end
0484
0485
0486 if plotflag
0487 if strcmp(method,'gpd,ml') | strcmp(method,'exp,ml')| strcmp(method,'exp,mld')
0488 Csum = cumsum(dN(:,2));
0489 q1 = ([1; Csum(1:end-1)+1; Csum(1:end)]-1/2)/n;
0490 xq1 = wgpdinv(q1,Est.k,Est.s)';
0491 x = [dN(:,1); dN(:,1)]+0.5;
0492 wqqplot(x,xq1); hold on
0493 plot([0 max(x)],[0 max(x)]), hold off
0494 xlabel('data'), ylabel('Quantile')
0495 elseif strcmp(methodShape,'gpd')
0496 wqqplot(x,xq1); hold on
0497 plot([0 max(x)],[0 max(x)]), hold off
0498 xlabel('data'), ylabel('Quantile')
0499 elseif strcmp(method,'ray,ml')
0500 Csum = cumsum(dN(:,2));
0501 q1 = ([1; Csum(1:end-1)+1; Csum(1:end)]-1/2)/n;
0502 xq1 = sqrt(offset^2-Est.s^2*log(1-q1))-offset;
0503 x = [dN(:,1); dN(:,1)]+0.5;
0504 wqqplot(x,xq1); hold on
0505 plot([0 max(x)],[0 max(x)]), hold off
0506 xlabel('data'), ylabel('Quantile')
0507 else
0508 plot(X(:,1),Y,'+'), hold on
0509 plot(X(:,1),X*th,'r'), hold off
0510 xlabel('data'), ylabel('log(lc)')
0511 end
0512 end
0513
0514
0515
0516