Home > wafo > trgauss > lc2tr2.m

lc2tr2

PURPOSE ^

Estimate transformation, g, from observed crossing intensity, version2.

SYNOPSIS ^

[g, test, g2] = lc2tr(cross,ma,sa,varargin);

DESCRIPTION ^

 LC2TR2 Estimate transformation, g, from observed crossing intensity, version2.
 
         Assumption: a Gaussian process, Y, is related to the
                     non-Gaussian process, X, by Y = g(X). 
 
   CALL:  [g, test,g2] = lc2tr2(lc,ma,sa,options);
 
      g,g2  = smoothed and empirical estimate of the transformation  g.     
      test  = test observator int (g(u)-u)^2 du  where int limits is
              given by param. This is a measure of departure of the 
              data from the Gaussian model.
      lc    = a two column matrix with levels in the first column
              and number of upcrossings in the second.
      ma,sa = mean and standard deviation of the process
 
    options = structure with the fields:
   csm,gsm  - defines the smoothing of the crossing intensity and the
              transformation g. Valid values must be
              0<=csm,gsm<=1. (default csm = 0.9 gsm=0.05)
              Smaller values gives smoother functions.
      param - vector which defines the region of variation of the data X.
              (default [-5 5 513]). 
   plotflag - 0 no plotting (Default)
              1 plots empirical and smoothed g(u) and theoretical for a Gaussian model. 
              2 monitor development of estimation
  linextrap - 0 use a regular smoothing spline 
              1 use a smoothing spline with a constraint on the ends to 
                ensure linear extrapolation outside the range of the data.
                (default)
  cvar      - Variances for the crossing intensity. (default 1)
  gvar      - Variances for the empirical transformation, g. (default  1) 
  ne        - Number of extremes (maxima & minima) to remove from the
               estimation of the transformation. This makes the
               estimation more robust against outliers. (default 7)
  Ntr        - Maximum length of empirical crossing intensity.
               The empirical crossing intensity is interpolated
               linearly  before smoothing if the length exceeds Ntr.
               A reasonable NTR will significantly speed up the
               estimation for long time series without loosing any
               accuracy. NTR should be chosen greater than
               PARAM(3). (default 1000)
 
     The empirical crossing intensity is usually very irregular.
   More than one local maximum of the empirical crossing intensity
   may cause poor fit of the transformation. In such case one
   should use a smaller value of GSM or set a larger variance for GVAR. 
     If X(t) is likely to cross levels higher than 5 standard deviations  
   then the vector param has to be modified.  For example if X(t) is 
   unlikely to cross a level of 7 standard deviations one can use 
   param = [-7 7 513].
 
  Example:
  Hm0 = 7;
  S = jonswap([],Hm0); g=ochitr([],[Hm0/4]); 
  S.tr=g;S.tr(:,2)=g(:,2)*Hm0/4;
  xs = spec2sdat(S,2^13);
  lc = dat2lc(xs);
  g0 = lc2tr2(lc,0,Hm0/4,'plot','iter');         % Monitor the development
  g1 = lc2tr2(lc,0,Hm0/4,troptset('gvar', .5 )); % Equal weight on all points
  g2 = lc2tr2(lc,0,Hm0/4,'gvar', [3.5 .5 3.5]);  % Less weight on the ends
  hold on, trplot(g1,g)                          % Check the fit
  trplot(g2)
 
  See also  troptset, dat2tr, trplot, findcross, smooth

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [g, test, g2] = lc2tr(cross,ma,sa,varargin);
002 %LC2TR2 Estimate transformation, g, from observed crossing intensity, version2.
003 %
004 %        Assumption: a Gaussian process, Y, is related to the
005 %                    non-Gaussian process, X, by Y = g(X). 
006 %
007 %  CALL:  [g, test,g2] = lc2tr2(lc,ma,sa,options);
008 %
009 %     g,g2  = smoothed and empirical estimate of the transformation  g.     
010 %     test  = test observator int (g(u)-u)^2 du  where int limits is
011 %             given by param. This is a measure of departure of the 
012 %             data from the Gaussian model.
013 %     lc    = a two column matrix with levels in the first column
014 %             and number of upcrossings in the second.
015 %     ma,sa = mean and standard deviation of the process
016 %
017 %   options = structure with the fields:
018 %  csm,gsm  - defines the smoothing of the crossing intensity and the
019 %             transformation g. Valid values must be
020 %             0<=csm,gsm<=1. (default csm = 0.9 gsm=0.05)
021 %             Smaller values gives smoother functions.
022 %     param - vector which defines the region of variation of the data X.
023 %             (default [-5 5 513]). 
024 %  plotflag - 0 no plotting (Default)
025 %             1 plots empirical and smoothed g(u) and theoretical for a Gaussian model. 
026 %             2 monitor development of estimation
027 % linextrap - 0 use a regular smoothing spline 
028 %             1 use a smoothing spline with a constraint on the ends to 
029 %               ensure linear extrapolation outside the range of the data.
030 %               (default)
031 % cvar      - Variances for the crossing intensity. (default 1)
032 % gvar      - Variances for the empirical transformation, g. (default  1) 
033 % ne        - Number of extremes (maxima & minima) to remove from the
034 %              estimation of the transformation. This makes the
035 %              estimation more robust against outliers. (default 7)
036 % Ntr        - Maximum length of empirical crossing intensity.
037 %              The empirical crossing intensity is interpolated
038 %              linearly  before smoothing if the length exceeds Ntr.
039 %              A reasonable NTR will significantly speed up the
040 %              estimation for long time series without loosing any
041 %              accuracy. NTR should be chosen greater than
042 %              PARAM(3). (default 1000)
043 %
044 %    The empirical crossing intensity is usually very irregular.
045 %  More than one local maximum of the empirical crossing intensity
046 %  may cause poor fit of the transformation. In such case one
047 %  should use a smaller value of GSM or set a larger variance for GVAR. 
048 %    If X(t) is likely to cross levels higher than 5 standard deviations  
049 %  then the vector param has to be modified.  For example if X(t) is 
050 %  unlikely to cross a level of 7 standard deviations one can use 
051 %  param = [-7 7 513].
052 %
053 % Example:
054 % Hm0 = 7;
055 % S = jonswap([],Hm0); g=ochitr([],[Hm0/4]); 
056 % S.tr=g;S.tr(:,2)=g(:,2)*Hm0/4;
057 % xs = spec2sdat(S,2^13);
058 % lc = dat2lc(xs);
059 % g0 = lc2tr2(lc,0,Hm0/4,'plot','iter');         % Monitor the development
060 % g1 = lc2tr2(lc,0,Hm0/4,troptset('gvar', .5 )); % Equal weight on all points
061 % g2 = lc2tr2(lc,0,Hm0/4,'gvar', [3.5 .5 3.5]);  % Less weight on the ends
062 % hold on, trplot(g1,g)                          % Check the fit
063 % trplot(g2)
064 %
065 % See also  troptset, dat2tr, trplot, findcross, smooth
066 
067 % NB! the transformated data will be N(0,1)
068 
069 % Reference
070 % Rychlik , I., Johannesson, P., and Leadbetter, M.R. (1997)
071 % "Modelling and statistical analysis of ocean wavedata 
072 % using a transformed Gaussian process",
073 % Marine structures, Design, Construction and Safety, 
074 % Vol 10, pp 13--47
075 % 
076 
077 
078 % Tested on: Matlab 5.3, 5.2, 5.1
079 % History:
080 % by pab 29.12.2000
081 % based on lc2tr, but the inversion is faster.
082 % by IR and PJ
083 
084 
085 opt = troptset('chkder','on','plotflag','off','csm',.9,'gsm',.05,....
086     'param',[-5 5 513],'delay',2,'linextrap','on','ntr',1000,'ne',7,'cvar',1,'gvar',1);
087 % If just 'defaults' passed in, return the default options in g
088 if nargin==1 & nargout <= 1 & isequal(cross,'defaults')
089   g = opt; 
090   return
091 end
092 error(nargchk(3,inf,nargin)) 
093 if nargin>=4 ,  opt  = troptset(opt,varargin{:}); end
094 csm2 = opt.gsm;
095 param = opt.param;
096 ptime = opt.delay;
097 Ne  = opt.ne;
098 switch opt.chkder;
099   case 'off', chkder = 0;
100   case 'on',  chkder = 1;
101   otherwise,  chkder = opt.chkder;
102 end
103 switch opt.linextrap;
104   case 'off', def = 0;
105   case 'on',  def = 1;
106   otherwise,  def = opt.linextrap;
107 end
108 switch opt.plotflag
109   case {'none','off'},   plotflag = 0;
110   case 'final', plotflag = 1;
111   case 'iter',  plotflag = 2;
112   otherwise,    plotflag = opt.plotflag;
113 end
114 ncr = length(cross);
115 if ncr>opt.ntr & opt.ntr>0,
116    x0 = linspace(cross(1+Ne,1),cross(end-Ne,1),opt.ntr)';
117    cros = [ x0,interp1q(cross(:,1),cross(:,2),x0)];
118    Ne = 0;
119    Ner = opt.ne;
120    ncr = opt.ntr;
121  else
122    Ner = 0;
123   cros=cross;
124 end
125 
126 ng = length(opt.gvar);
127 if ng==1
128   gvar = opt.gvar(ones(ncr,1));
129 else
130   gvar = interp1(linspace(0,1,ng)',opt.gvar(:),linspace(0,1,ncr)','*linear');  
131 end
132 ng = length(opt.cvar);
133 if ng==1
134   cvar = opt.cvar(ones(ncr,1));
135 else
136   cvar = interp1(linspace(0,1,ng)',opt.cvar(:),linspace(0,1,ncr)','*linear');  
137 end
138 
139 g = zeros(param(3),2);
140 
141 uu = levels(param);
142 
143 g(:,1) = sa*uu' + ma;
144 
145 g2   = cros;
146 
147 if Ner>0, % Compute correction factors
148  cor1 = trapz(cross(1:Ner+1,1),cross(1:Ner+1,2));
149  cor2 = trapz(cross(end-Ner-1:end,1),cross(end-Ner-1:end,2));
150 else
151   cor1 = 0;
152   cor2 = 0;
153 end
154 cros(:,2) = cumtrapz(cros(:,1),cros(:,2))+cor1;
155 cros(:,2) = (cros(:,2)+.5)/(cros(end,2) + cor2 +1);
156 cros(:,1) = (cros(:,1)-ma)/sa;
157 
158 % find the mode
159 [tmp,imin]= min(abs(cros(:,2)-.15));
160 [tmp,imax]= min(abs(cros(:,2)-.85));
161 inde = imin:imax;
162 tmp =  smooth(cros(inde,1),g2(inde,2),opt.csm,cros(inde,1),def,cvar(inde));
163 
164 [tmp imax] = max(tmp);
165 u0 = cros(inde(imax),1);
166 %u0 = interp1q(cros(:,2),cros(:,1),.5)
167 
168 
169 cros(:,2) = wnorminv(cros(:,2),-u0,1);
170 
171 g2(:,2)   = cros(:,2);
172 % NB! the smooth function does not always extrapolate well outside the edges
173 % causing poor estimate of g  
174 % We may alleviate this problem by: forcing the extrapolation
175 % to be linear outside the edges or choosing a lower value for csm2.
176 
177 inds = 1+Ne:ncr-Ne;% indices to points we are smoothing over
178 scros2 = smooth(cros(inds,1),cros(inds,2),csm2,uu,def,gvar(inds));
179 
180 g(:,2) = scros2';%*sa; %multiply with stdev 
181 
182 if chkder~=0
183    for ix = 1:5
184     dy = diff(g(:,2));
185     if any(dy<=0)
186       warning('The empirical crossing spectrum is not sufficiently smoothed.')
187       disp('        The estimated transfer function, g, is not ')
188       disp('        a strictly increasing function.')
189       dy(dy>0)=eps;
190       gvar = -([dy;0]+[0;dy])/2+eps;
191       g(:,2) = smooth(g(:,1),g(:,2),1,g(:,1),def,ix*gvar);
192     else 
193       break
194     end
195   end
196 end
197 if 0,  %either
198   test = sqrt((param(2)-param(1))/(param(3)-1)*sum((uu-scros2).^2));
199 else % or
200   %test=sqrt(simpson(uu,(uu-scros2).^2));
201 % or
202   test=sqrt(trapz(uu,(uu-scros2).^2));
203 end 
204 
205 
206 if plotflag>0, 
207   trplot(g ,g2,ma,sa)
208   %legend(['Smoothed (T='  num2str(test) ')'],'g(u)=u','Not smoothed',0)
209   %ylabel('Transfer function g(u)')
210   %xlabel('Crossing level u')
211   
212   if plotflag>1,pause(ptime),end
213 end
214 
215 
216 
217 
218 
219 
220 
221

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