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: 
    x = load('sea.dat'); 
    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*') 
  
  See also  rfmextrapoalte

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

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: 
039 %   x = load('sea.dat'); 
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 % 
046 % See also  rfmextrapoalte 
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. 
057 %   Added GPD. 
058 % Updated by PJ  29-Oct-2003 
059 %   If N<1, estimate parameters Pout, but set tpe=[]. 
060 % Updated by PJ  11-Mar-2005 
061 %   Added "Conservative" extrapolation of load spectrum. 
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