Return period value with POT method using WAFO
ke
Threshold method. I am seeking the expected max value of a structural
load during a 3-hour storm, based on a 3-hour time series (in real
world units) from a wave tank. The load reaches a peak as each wave
passes by. I am attempting to follow the method in
www.marine.csiro.au/~hem021/reports/extremes.pdf . Is this a sensible
approach,
1) In Matlab, identify peak values after declustering. Choose
threshold load value ('Threshold'). Find excess of the peaks over
threshold ('PeakExcess').
2) In WAFO, fit GPD distribution and check that it matches observed
distribution,
gpd = wgpdfit(PeakExcess, 'pwm', 1) ;
k_hat = gpd(1) ; % GPD scale parameter
alpha_hat = gpd(2) ; % GPD shape parameter
3) Verify fit with quantile-quantile plot,
% Generate random variable with GPD distribution
R = wgpdrnd(k_hat, alpha_hat, 0, [ length(PeakExcess) 1]) ;
% Compare to observed distribution
wqqplot(PeakExcess, R) ; % Should be a straight line
4) Estimate number of clusters/peaks per return period
ReturnPeriodSeconds = 3600*3 ; % Duration of return period in
seconds
TestDuration = range(t) ; % Duration of the wave tank test in real-
world seconds
NumberPeaksPerReturnPeriod =
ReturnPeriodSeconds*length(PeakExcess)./TestDuration ;
5) Estimate return period value for GPD
if k_hat ~= 0
ValueAtReturnPeriod = Threshold + (alpha_hat/k_hat)*(1 -
NumberPeaksPerReturnPeriod^-k_hat) ;
end
Thanks for any thoughts or warnings!
Per Andreas Brodtkorb
http://code.google.com/p/wafo/source/browse/trunk/wafo25/statistics/findpot.m
http://code.google.com/p/wafo/source/browse/trunk/wafo25/statistics/decluster.m
See also:
http://code.google.com/p/wafo/source/browse/trunk/wafo25/statistics/extremalidx.m
http://code.google.com/p/wafo/source/browse/trunk/wafo25/statistics/extremalidxrange.m
In order to select a good threshold you may also use:
http://code.google.com/p/wafo/source/browse/trunk/wafo25/statistics/disprsnidx.m
http://code.google.com/p/wafo/source/browse/trunk/wafo25/statistics/disprsnidxrange.m
With WAFO 2.5 you could do something like this:
Nmin=7; % minimum number of extremes
Tmin=5; % minimum distance between extremes
Tb = 15; % block period
xn = load('sea.dat');
timeSpan = (xn(end,1)-xn(1,1))/60; % in minutes
dt = xn(2,1)-xn(1,1); % in seconds
tc = dat2tc(xn);
umin = median(tc(tc(:,2)>0,2));
Ie0 = findpot(tc, 0.9*umin, Tmin);
Ev = sort(tc(Ie0,2));
Ne = numel(Ev)
if Ne>Nmin && Ev(Ne-Nmin)>umin,
umax = Ev(Ne-Nmin);
else
umax = umin;
end
Nu = floor((umax-umin)/0.025)+1;
u = linspace(umin,umax,Nu);
mrl = reslife(Ev, 'u',u);
umin0 = umin;
for io= numel(mrl.data):-1:1,
CI = mrl.dataCI(io:end,:);
if ~(max(CI(:,1))<= mrl.data(io) & mrl.data(io)<=min(CI(:,2))),
umin0 = mrl.args(io);
break;
end
end
[di, threshold, ok_u] = disprsnidx(tc(Ie0,:), 'Tb', Tb, 'alpha',0.05, 'u',u);
figure(1)
plot(di)
vline(threshold) % threshold from dispersion index
vline(umin0,'g') % compared to the one obtained from mean residual life plot.
figure(2)
plot(mrl)
vline(threshold) % threshold from dispersion index
vline(umin0,'g') % compared to the one obtained from mean residual life plot.
% threshold around 1.2 seems appropriate
Ie = findpot(tc, threshold, Tmin);
lambda = numel(Ie)/timeSpan; % # Y> threshold per minute
varLambda = lambda*(1-(dt/60)*lambda)/timeSpan;
stdLambd = sqrt(varLambda)
Ev = tc(Ie,2);
phat = fitgenpar(Ev, 'fixpar',[nan,nan,threshold], method,'mps');
figure(3)
phat.plotfitsumry() % check fit to data
% Return period and return level
Tr = 60 %Return period in minutes
[xr,xrlo,xrup] = invgenpar(1./(lambda*Tr),phat,'lowertail',false,
'alpha', 0.05) % 60 minutes return level + 95%CI
[xr,xrlo5,xrup5] = invgenpar(1./(lambda*Tr),phat,'lowertail',false,
'alpha', 0.5) % 60 minutes return level + 50%CI
Tr = 3*60 %Return period in minutes
[xr,xrlo,xrup] = invgenpar(1./(lambda*Tr),phat,'lowertail',false,
'alpha', 0.05) % 180 minutes return level + 95%CI
[xr,xrlo5,xrup5] = invgenpar(1./(lambda*Tr),phat,'lowertail',false,
'alpha', 0.5) % 180 minutes return level + 50%CI
Best regards Per A.
> --
> Du mottar denne meldingen fordi du abonnerer på Google-gruppen «wafo».
> Hvis du vil legge inn en melding i denne gruppen, kan du sende e-post til wa...@googlegroups.com.
> Hvis du vil avslutte abonnementet på denne gruppen, sender du en e-post til wafo+uns...@googlegroups.com.
> Hvis du vil ha flere alternativer, kan du besøke gruppen på http://groups.google.com/group/wafo?hl=no.
>
>
ke
On Feb 28, 5:39 am, Per Andreas Brodtkorb
<per.andreas.brodtk...@gmail.com> wrote:
> In order to decluster the data you may use:http://code.google.com/p/wafo/source/browse/trunk/wafo25/statistics/f...http://code.google.com/p/wafo/source/browse/trunk/wafo25/statistics/d...
>
> See also:http://code.google.com/p/wafo/source/browse/trunk/wafo25/statistics/e...http://code.google.com/p/wafo/source/browse/trunk/wafo25/statistics/e...
>
> In order to select a good threshold you may also use:http://code.google.com/p/wafo/source/browse/trunk/wafo25/statistics/d...http://code.google.com/p/wafo/source/browse/trunk/wafo25/statistics/d...
>
> > I would like to calculate a return period value using the Peak Over
> > Threshold method. I am seeking the expected max value of a structural
> > load during a 3-hour storm, based on a 3-hour time series (in real
> > world units) from a wave tank. The load reaches a peak as each wave
> > passes by. I am attempting to follow the method in
ke
calculation (disprsnidx) is Tb=15 seconds? There are only about 3
peaks within each 15-second interval, which does not seem like enough
to test whether they are distributed as a Poisson process. But I am
sure I am misunderstanding Tb.
On Feb 28, 5:39 am, Per Andreas Brodtkorb
<per.andreas.brodtk...@gmail.com> wrote:
> In order to decluster the data you may use:http://code.google.com/p/wafo/source/browse/trunk/wafo25/statistics/f...http://code.google.com/p/wafo/source/browse/trunk/wafo25/statistics/d...
>
> See also:http://code.google.com/p/wafo/source/browse/trunk/wafo25/statistics/e...http://code.google.com/p/wafo/source/browse/trunk/wafo25/statistics/e...
>
>
> > I would like to calculate a return period value using the Peak Over
> > Threshold method. I am seeking the expected max value of a structural
> > load during a 3-hour storm, based on a 3-hour time series (in real
> > world units) from a wave tank. The load reaches a peak as each wave
> > passes by. I am attempting to follow the method in
Per Andreas Brodtkorb
ke
my time series are independent, correct?
On Mar 16, 6:50 pm, Per Andreas Brodtkorb
> The reason for selecting Tb=15s is that the autocorrelation for the process
> is not significantly different from zero in this example. Assuming that the
> autocorrelation is zero for t>15s and that you consider the data as gaussian
> distributed, it means that events further than 15s apart are independent.
> However, the assumptions might not be correct so that another choice may be
> more suitable.
>
>
> > Can you please explain why the block period for the Dispersion Index
> > calculation (disprsnidx) is Tb=15 seconds? There are only about 3
> > peaks within each 15-second interval, which does not seem like enough
> > to test whether they are distributed as a Poisson process. But I am
> > sure I am misunderstanding Tb.
>
> > On Feb 28, 5:39 am, Per Andreas Brodtkorb
> > <per.andreas.brodtk...@gmail.com> wrote:
> > > In order to decluster the data you may use:
> > ..
>
> > > See also:
> >http://code.google.com/p/wafo/source/browse/trunk/wafo25/statistics/e....
> > ..
> > > In order to select a good threshold you may also use:
> > ..