Home > wafo > cycles > tpextrapolate.m

# tpextrapolate

## PURPOSE

Extrapolates a sequence of turning points.

## SYNOPSIS

[tpe,Pout,I,tpe0] = tpextrapolate(tp,N,Pin,plotflag)

## DESCRIPTION

``` TPEXTRAPOLATE Extrapolates a sequence of turning points.

CALL: tpe = tpextrapolate(tp,N);
[tpe,Pout,I] = tpextrapolate(tp,N,Pin,plotflag)

tp          = A sequence of turning points.         [n,1] / [n,2]
N           = Number of blocks to extrapolate.
Pin         = Input parameters. (optional)          [struct array]
.method    = Method for extrapolating.
'exp' : Exponential distribution  (default)
'gpd' : Generalized Pareto Distribution.
Note: Use 'gpd' only if 'exp' gives bad fit to data.
.LCfrac    = Fraction of level crossings (LC) for automatic choice of levels, u_lev.
The upper and lower levels for extrapolation are chosen so that they are
crossed LCmax*LCfrac number of times.  LCmax is the maximum of LC.
Default: Pin.LCfrac=1/sqrt(LCmax)
.u_lev     = Lower and upper levels for extrapolation.  (manual choice of levels)
Pin.u_lev=[u_min u_max]
.lim       = Limits of signal,  Pin.lim=[minlim maxlim] .
Values outside the range is set to the limit.
Default: Pin.lim=[], no limits

plotflag    = 0: Don't plot diagnostic plots, (default)
1: Plot exponential probability plots,
2: Also plot mean excess plots.

tpe         = Extrapolated turning points           [N*n,1] / [N*n,2]
Pout        = Output parameters. (see Pin)                   [struct array]
I           = Indeces to the extrapolated points.

The highest maxima and the lowest minima of the turning points are extrapolated.
The result is an N times as long signal consisting of N blocks.  Each block is
generated from tp with randomly simulated maxima above u_max, and minima below u_min.
The mean of the exceedances above u_max and below u_min are estimated from tp.

Example:
tp = dat2tp(x,0.5);
tpe = tpextrapolate(tp,1,[],1);
clf, plot(tp(:,1),tp(:,2),'b',tpe(:,1),tpe(:,2),'r')
[tpe,Pout,I] = tpextrapolate(tp,1,[],2);
clf, plot(tpe(:,1),tpe(:,2),'b',tpe(I.min,1),tpe(I.min,2),'g*',tpe(I.max,1),tpe(I.max,2),'g*')

## CROSS-REFERENCE INFORMATION

This function calls:
 rfcfilter Rainflow filter a signal. tp2lc Calculates the number of upcrossings from the turning points. wexpplot Plots data on a Exponential distribution paper wexprnd Random matrices from an Exponential distribution wgpdfit Parameter estimates for Generalized Pareto data wgpdrnd Random matrices from a Generalized Pareto Distribution drawnow Flush pending graphics events. error Display message and abort function. fieldnames Get structure field names. figure Create figure window. getfield Get structure field contents. grid Grid lines. isnumeric True for numeric arrays. linspace Linearly spaced vector. lower Convert string to lowercase. mean Average or mean value. num2str Convert number to string. (Fast version) plot Linear plot. setfield Set structure field contents. strcmp Compare strings. subplot Create axes in tiled positions. title Graph title. xlabel X-axis label. ylabel Y-axis label.
This function is called by:
 test_cycles Quick test of the routines in module 'cycles'

## SOURCE CODE

```001 function [tpe,Pout,I,tpe0] = tpextrapolate(tp,N,Pin,plotflag)
002
003 %TPEXTRAPOLATE Extrapolates a sequence of turning points.
004 %
005 % CALL: tpe = tpextrapolate(tp,N);
006 %       [tpe,Pout,I] = tpextrapolate(tp,N,Pin,plotflag)
007 %
008 %   tp          = A sequence of turning points.         [n,1] / [n,2]
009 %   N           = Number of blocks to extrapolate.
010 %   Pin         = Input parameters. (optional)          [struct array]
011 %    .method    = Method for extrapolating.
012 %                'exp' : Exponential distribution  (default)
013 %                'gpd' : Generalized Pareto Distribution.
014 %                        Note: Use 'gpd' only if 'exp' gives bad fit to data.
015 %    .LCfrac    = Fraction of level crossings (LC) for automatic choice of levels, u_lev.
016 %                 The upper and lower levels for extrapolation are chosen so that they are
017 %                 crossed LCmax*LCfrac number of times.  LCmax is the maximum of LC.
018 %                 Default: Pin.LCfrac=1/sqrt(LCmax)
019 %    .u_lev     = Lower and upper levels for extrapolation.  (manual choice of levels)
020 %                   Pin.u_lev=[u_min u_max]
021 %    .lim       = Limits of signal,  Pin.lim=[minlim maxlim] .
022 %                 Values outside the range is set to the limit.
023 %                 Default: Pin.lim=[], no limits
024 %
025 %   plotflag    = 0: Don't plot diagnostic plots, (default)
026 %                 1: Plot exponential probability plots,
027 %                 2: Also plot mean excess plots.
028 %
029 %   tpe         = Extrapolated turning points           [N*n,1] / [N*n,2]
030 %   Pout        = Output parameters. (see Pin)                   [struct array]
031 %   I           = Indeces to the extrapolated points.
032 %
033 % The highest maxima and the lowest minima of the turning points are extrapolated.
034 % The result is an N times as long signal consisting of N blocks.  Each block is
035 % generated from tp with randomly simulated maxima above u_max, and minima below u_min.
036 % The mean of the exceedances above u_max and below u_min are estimated from tp.
037 %
038 % Example:
040 %   tp = dat2tp(x,0.5);
041 %   tpe = tpextrapolate(tp,1,[],1);
042 %   clf, plot(tp(:,1),tp(:,2),'b',tpe(:,1),tpe(:,2),'r')
043 %   [tpe,Pout,I] = tpextrapolate(tp,1,[],2);
044 %   clf, plot(tpe(:,1),tpe(:,2),'b',tpe(I.min,1),tpe(I.min,2),'g*',tpe(I.max,1),tpe(I.max,2),'g*')
045 %
047
048 % Tested  on Matlab  6.5
049 %
050 % History:
051 % Created by PJ (Pär Johannesson) 16-Apr-2003
052 % Updated by PJ  11-Jun-2003
053 % Updated by PJ  24-Jun-2003
054 % Updated by PJ  05-Sep-2003
055 %   Added output I, and input Pin.lim.
056 %   Now also handles zero number of exceedances.
058 % Updated by PJ  29-Oct-2003
059 %   If N<1, estimate parameters Pout, but set tpe=[].
060 % Updated by PJ  11-Mar-2005
062
063
064 % Check input arguments
065 ni = nargin;
066 no = nargout;
067 error(nargchk(2,4,ni));
068
069 if ni<3, Pin=[]; end
070 if ni<4, plotflag=[]; end
071
072 if isempty(plotflag)
073     plotflag=0;
074 end
075
076 %
077 % Default Parameters
078 %
079
080 Pout.method = 'exp';
081 Pout.LCfrac = []; %0.05;
082 Pout.u_lev = [];
083 Pout.lim = [];
084
085 % Copy input parameters Pin to Pout
086 if ~isempty(Pin)
087     Fname = fieldnames(Pin);
088     for i = 1:length(Fname)
089         Pout = setfield(Pout,Fname{i},getfield(Pin,Fname{i}));
090     end
091 end
092
093 method = lower(Pout.method);
094 if ~(strcmp(method,'exp') | strcmp(method,'gpd'))
095     error(['Undefined method: ' Pout.method]),
096 end
097
098 if isempty(Pout.u_lev)
099     lc = tp2lc(tp);  % Level crossings
100
101     LCmax=max(lc(:,2));
102     if isempty(Pout.LCfrac)
103         Pout.LCfrac = 1/sqrt(LCmax);
104     end
105
106     nLCfrac = LCmax*Pout.LCfrac;
107     imax = max(find(lc(:,2)>nLCfrac));
108     imin = min(find(lc(:,2)>nLCfrac));
109     umax = lc(imax,1);
110     umin = lc(imin,1);
111     Pout.u_lev = [umin umax];
112 end
113
114 u_min = Pout.u_lev(1);
115 u_max = Pout.u_lev(2);
116
117 if tp(1,end)>tp(2,end)
118     StartMax=1; StartMin=2;
119 else
120     StartMax=2; StartMin=1;
121 end
122
123 % Estimate excedances
124
125 Imax = 2*(find(tp(StartMax:2:end,end)>u_max)-1)+StartMax;
126 Imin = 2*(find(tp(StartMin:2:end,end)<u_min)-1)+StartMin;
127
128 y_max = tp(Imax,end)-u_max;
129 y_min = -(tp(Imin,end)-u_min);
130
131 a_max = mean(y_max);
132 a_min = mean(y_min);
133 if strcmp(method,'gpd')
134     [gpd_max,cov_max] = wgpdfit(y_max,'ml',0);
135     [gpd_min,cov_min] = wgpdfit(y_min,'ml',0);
136     ci_k_max = [gpd_max(1)-2*sqrt(cov_max(1)) gpd_max(1)+2*sqrt(cov_max(1))]
137     ci_k_min = [gpd_min(1)-2*sqrt(cov_min(1)) gpd_min(1)+2*sqrt(cov_min(1))]
138 end
139
140 if plotflag>0
141     subplot(2,2,1),plot(1:length(y_max),y_max,'.'),
142     title(['Exceedances of maxima above u_{max}=' num2str(u_max)]), ylabel('Exceedances of maxima')
143     subplot(2,2,2), if length(y_max)>0, wexpplot(y_max), end
144     subplot(2,2,3),plot(1:length(y_min),y_min,'.'),
145     title(['Exceedances of minima below u_{min}=' num2str(u_min)]), ylabel('Exceedances of minima')
146     subplot(2,2,4), if length(y_min)>0, wexpplot(y_min), end
147     drawnow
148
149     Pout.Zmin = y_min;
150     Pout.Zmax = y_max;
151 end
152
153 if N>0
154     meth=2;  % Method 2
155     conservarive = 1;  % 1=Conservative extrapolation
156
157     if meth==2
158         [yy_max,IImax] = sort(y_max);
159         [yy_min,IImin] = sort(y_min);
160     end
161     tpe = [];
162     for k = 1:N
163
164         % Simulate independent Exp or GPD exceedances
165
166         if length(y_max)>0
167             if strcmp(method,'exp')
168                 yr_max = wexprnd(a_max,length(y_max),1);
169             else
170                 yr_max = wgpdrnd(gpd_max(1),gpd_max(2),0,length(y_max),1);
171             end
172         else
173             yr_max=[];
174         end
175
176         if length(y_min)>0
177             if strcmp(method,'exp')
178                 yr_min = wexprnd(a_min,length(y_min),1);
179             else
180                 yr_min = wgpdrnd(gpd_min(1),gpd_min(2),0,length(y_min),1);
181             end
182         else
183             yr_min=[];
184         end
185
186         if meth ==1
187             % Method 1
188             % Independent ordering of excedances
189             tpr = tp;
190             tpr(Imax,end) = u_max + yr_max;
191             tpr(Imin,end) = u_min - yr_min;
192         else
193
194             % Method 2
195             % Simulate independent Exp
196             % Order the sample according to the order of the measurements.
197
198             yr_max = sort(yr_max);
199             yr_min = sort(yr_min);
200
201             tpr = tp;
202             tpr(Imax(IImax),end) = u_max+yr_max;
203             tpr(Imin(IImin),end) = u_min-yr_min;
204         end
205
206         if conservarive
207             I=find(tpr(Imin,end)>tp(Imin,end));
208             tpr(Imin(I),end) = tp(Imin(I),end);
209             I=find(tpr(Imax,end)<tp(Imax,end));
210             tpr(Imax(I),end) = tp(Imax(I),end);
211         end
212
213         tpe = [tpe; tpr];
214     end
215
216     if no>3
217         tpe0=tpe;
218     end
219
220     % Apply limits
221     if ~isempty(Pout.lim)
222         minlim=Pout.lim(1);  maxlim=Pout.lim(2);
223         if isnumeric(maxlim) % Don't apply if maxlim=Inf
224             I = find(tpe(:,end)>maxlim);
225             if ~isempty(I)
226                 tpe(I,end)=maxlim;
227             end
228         end
229         if isnumeric(minlim) % Don't apply if minlim=-Inf
230             I = find(tpe(:,end)<minlim);
231             if ~isempty(I)
232                 tpe(I,end)=minlim;
233             end
234         end
235     end
236
237     % Make sure that the output is a sequence of turning points
238     tpe = rfcfilter(tpe,0,1);
239     %tpe = dat2tp(tpe);
240
241 else
242     tpe = [];
243 end
244
245 % Output parameters
246 Pout.a_min = a_min;
247 Pout.a_max = a_max;
248 Pout.n_min = length(y_min);
249 Pout.n_max = length(y_max);
250
251 if strcmp(method,'gpd')
252     Pout.gpd_min = gpd_min;
253     Pout.gpd_max = gpd_max;
254     Pout.ci_k_min = ci_k_min;
255     Pout.ci_k_max = ci_k_max;
256 end
257
258 if N>0
259     % Indeces of extrapolated values
260     if no>2
261         if tpe(1,end)>tpe(2,end)
262             StartMax=1; StartMin=2;
263         else
264             StartMax=2; StartMin=1;
265         end
266
267         % Estimate excedances
268         I=[];
269         I.max = 2*(find(tpe(StartMax:2:end,end)>u_max)-1)+StartMax;
270         I.min = 2*(find(tpe(StartMin:2:end,end)<u_min)-1)+StartMin;
271     end
272 end
273
274 % Diagnostic plots 2
275 if plotflag>1
276
277     m = mean(tp(:,end));
278
279     n=500;
280     Umin = linspace(min(tp(:,end)),m,n);
281     Umax = linspace(m,max(tp(:,end)),n);
282     Amin=zeros(1,n); Amax=zeros(1,n);
283     dAmin=zeros(1,n); dAmax=zeros(1,n);
284     for k=1:n
285         Imax = 2*(find(tp(StartMax:2:end,end)>Umax(k))-1)+StartMax;
286         Imin = 2*(find(tp(StartMin:2:end,end)<Umin(k))-1)+StartMin;
287
288         y_max = tp(Imax,end)-Umax(k);
289         if ~isempty(y_max), Amax(k) = mean(y_max); dAmax(k)=Amax(k)/sqrt(length(y_max)); end
290         y_min = -(tp(Imin,end)-Umin(k));
291         if ~isempty(y_min), Amin(k) = mean(y_min); dAmin(k)=Amin(k)/sqrt(length(y_min)); end
292     end
293
294     figure
295     subplot(2,1,1)
296     plot(Umax,Amax,'r',Umax,Amax-2*dAmax,'b:',Umax,Amax+2*dAmax,'b:'), grid
297     title('Maxima - Choose a level where the estimate is stable'),
298     ylabel('Estimated mean exceedance, m'), xlabel('Upper threshold level, u_{max}')
299     subplot(2,1,2)
300     plot(-Umin,Amin,'r',-Umin,Amin-2*dAmin,'b:',-Umin,Amin+2*dAmin,'b:'), grid
301     title('Minima - Choose a level where the estimate is stable'),
302     ylabel('Estimated mean exceedance, m'), xlabel('Lower threshold level, -u_{min}')
303
304     Pout.Umin = Umin;
305     Pout.Amin = Amin;
306     Pout.Umax = Umax;
307     Pout.Amax = Amax;
308 end
309
310```

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