Home > wafo > trgauss > lc2tr.m

lc2tr

PURPOSE ^

Estimate transformation, g, from observed crossing intensity.

SYNOPSIS ^

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

DESCRIPTION ^

 LC2TR Estimate transformation, g, from observed crossing intensity.
 
         Assumption: a Gaussian process, Y, is related to the
                     non-Gaussian process, X, by Y = g(X). 
        
   CALL:  [g, test,g2] = lc2tr(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 logarithm of crossing intensity 
              and the transformation g, respectively. 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 the theoretical for a Gaussian model. 
              2 monitor the development of the 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 logarithm of 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)
 
     The empirical crossing intensity is usually very irregular.
   More than one local maximum of the smoothed crossing intensity
   may cause poor fit of the transformation. In such case one
   should use a smaller value of CSM or set a larger variance for CVAR. 
     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 = lc2tr(lc,0,Hm0/4,'plot','iter');         % Monitor the development
  g1 = lc2tr(lc,0,Hm0/4,troptset('gvar', .5 )); % Equal weight on all points
  g2 = lc2tr(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 ^

0001 function [g, test, g2] = lc2tr(cross,ma,sa,varargin);
0002 %LC2TR Estimate transformation, g, from observed crossing intensity.
0003 %
0004 %        Assumption: a Gaussian process, Y, is related to the
0005 %                    non-Gaussian process, X, by Y = g(X). 
0006 %       
0007 %  CALL:  [g, test,g2] = lc2tr(lc,ma,sa,options);
0008 %
0009 %     g,g2  = smoothed and empirical estimate of the transformation  g.     
0010 %     test  = test observator int (g(u)-u)^2 du  where int limits is
0011 %             given by param. This is a measure of departure of the 
0012 %             data from the Gaussian model.
0013 %     lc    = a two column matrix with levels in the first column
0014 %             and number of upcrossings in the second.
0015 %     ma,sa = mean and standard deviation of the process
0016 %
0017 %   options = structure with the fields:
0018 % csm, gsm  - defines the smoothing of the logarithm of crossing intensity 
0019 %             and the transformation g, respectively. Valid values must 
0020 %             be 0<=csm,gsm<=1. (default csm=0.9, gsm=0.05)
0021 %             Smaller values gives smoother functions.
0022 %  param    - vector which defines the region of variation of the data X.
0023 %             (default [-5 5 513]). 
0024 %  plotflag - 0 no plotting (Default)
0025 %             1 plots empirical and smoothed g(u) and the theoretical for a Gaussian model. 
0026 %             2 monitor the development of the estimation
0027 % linextrap - 0 use a regular smoothing spline 
0028 %             1 use a smoothing spline with a constraint on the ends to 
0029 %               ensure linear extrapolation outside the range of the data.
0030 %               (default)
0031 % cvar      - Variances for the logarithm of the crossing intensity. (default  1) 
0032 % gvar      - Variances for the empirical transformation, g. (default  1) 
0033 % ne        - Number of extremes (maxima & minima) to remove from the
0034 %              estimation of the transformation. This makes the
0035 %              estimation more robust against outliers. (default 7)
0036 %
0037 %    The empirical crossing intensity is usually very irregular.
0038 %  More than one local maximum of the smoothed crossing intensity
0039 %  may cause poor fit of the transformation. In such case one
0040 %  should use a smaller value of CSM or set a larger variance for CVAR. 
0041 %    If X(t) is likely to cross levels higher than 5 standard deviations  
0042 %  then the vector param has to be modified.  For example if X(t) is 
0043 %  unlikely to cross a level of 7 standard deviations one can use 
0044 %  param = [-7 7 513].
0045 %
0046 % Example:
0047 % Hm0 = 7;
0048 % S = jonswap([],Hm0); g=ochitr([],[Hm0/4]); 
0049 % S.tr=g;S.tr(:,2)=g(:,2)*Hm0/4;
0050 % xs = spec2sdat(S,2^13);
0051 % lc = dat2lc(xs);
0052 % g0 = lc2tr(lc,0,Hm0/4,'plot','iter');         % Monitor the development
0053 % g1 = lc2tr(lc,0,Hm0/4,troptset('gvar', .5 )); % Equal weight on all points
0054 % g2 = lc2tr(lc,0,Hm0/4,'gvar', [3.5 .5 3.5]);  % Less weight on the ends
0055 % hold on, trplot(g1,g)                                   % Check the fit
0056 % trplot(g2)
0057 %
0058 % See also  troptset, dat2tr, trplot, findcross, smooth
0059 
0060 % NB! the transformated data will be N(0,1)
0061 
0062 % Reference
0063 % Rychlik , I., Johannesson, P., and Leadbetter, M.R. (1997)
0064 % "Modelling and statistical analysis of ocean wavedata 
0065 % using a transformed Gaussian process",
0066 % Marine structures, Design, Construction and Safety, 
0067 % Vol 10, pp 13--47
0068 % 
0069 
0070 
0071 % Tested on: Matlab 5.3, 5.2, 5.1
0072 % History:
0073 % revised pab 21.12.2000
0074 % - added interp of cross -> much faster estimation
0075 % - added example, chkder
0076 % - replaced optional arguments with a options struct
0077 % - added options to input: ne, cvar,gvar
0078 % - changed names: csm1 to csm
0079 %                  csm2 to gsm
0080 % - replaced monitor with plotflag==2
0081 % - default param is now [-5 5 513] -> better to have the discretization
0082 %  represented with exact numbers, especially when calculating
0083 %  derivatives of the transformation numerically.
0084 %
0085 %  modified by svi 29.09.99
0086 %  Transformations  g2, g2 are not normalized any longer.
0087 % revised pab 11.08.99
0088 % changed name from cross2tr to lc2tr
0089 %
0090 % modified by Per A. Brodtkorb 15.08.98
0091 %  to check if the smoothing is sufficient and 
0092 %  changed the calculation of the test statistic.
0093 %  moved the plotting routine to trplot 
0094 
0095 %  Beware of the problem that Carl de Boor's smooth function 
0096 %  does not always extrapolate well outside the ends when 
0097 %  the smoothing parameter, p, is close to one.
0098 %  Particularly extrapolation in the first smoothing may corrupt 
0099 %  the estimate of g. One solution to the problem is to 
0100 %  extrapolate linearly. This is incorporated into the smooth function.
0101 %  Yet, another solution is choosing a lower value for csm1
0102 %  or not to extrapolate at all in the first smoothing but instead leave
0103 %  all the extrapolation to the second smoothing. 
0104 %  (Probably better since csm2<<csm1))
0105 %
0106 %   also added secret options: plotflag and monitor the steps of
0107 %   estimation of the transformation 
0108 %
0109 
0110 opt = troptset('chkder','on','plotflag','off','csm',.95,'gsm',.05,....
0111     'param',[-5 5 513],'delay',2,'ntr',1000,'linextrap','on','ne',7,'cvar',1,'gvar',1);
0112 % If just 'defaults' passed in, return the default options in g
0113 if nargin==1 & nargout <= 1 & isequal(cross,'defaults')
0114   g = opt; 
0115   return
0116 end
0117 error(nargchk(3,inf,nargin)) 
0118 if nargin>=4,  opt  = troptset(opt,varargin{:}); end
0119 
0120 csm1 = opt.csm;
0121 csm2 = opt.gsm;
0122 param = opt.param;
0123 ptime = opt.delay;
0124 Ne  = opt.ne;
0125 switch opt.chkder;
0126   case 'off', chkder = 0;
0127   case 'on',  chkder = 1;
0128   otherwise,  chkder = opt.chkder;
0129 end
0130 switch opt.linextrap;
0131   case 'off', def = 0;
0132   case 'on',  def = 1;
0133   otherwise,  def = opt.linextrap;
0134 end
0135 switch opt.plotflag
0136   case {'none','off'},   plotflag = 0;
0137   case 'final', plotflag = 1;
0138   case 'iter',  plotflag = 2;
0139   otherwise,    plotflag = opt.plotflag;
0140 end
0141 ncr = length(cross);
0142 if ncr>opt.ntr & opt.ntr>0,
0143    x0 = linspace(cross(1+Ne,1),cross(end-Ne,1),opt.ntr)';
0144    cros = [ x0,interp1q(cross(:,1),cross(:,2),x0)];
0145    Ne=0;
0146    ncr = opt.ntr;
0147 else
0148   cros=cross;
0149 end
0150 
0151 
0152 ng = length(opt.gvar);
0153 if ng==1
0154   gvar = opt.gvar(ones(ncr,1));
0155 else
0156   gvar = interp1(linspace(0,1,ng)',opt.gvar(:),linspace(0,1,ncr)','*linear');  
0157 end
0158 ng = length(opt.cvar);
0159 if ng==1
0160   cvar = opt.cvar(ones(ncr,1));
0161 else
0162   cvar = interp1(linspace(0,1,ng)',opt.cvar(:),linspace(0,1,ncr)','*linear');  
0163 end  
0164 
0165 
0166 g = zeros(param(3),2);
0167 
0168 uu = levels(param);
0169 
0170 g(:,1) = sa*uu' + ma;
0171 
0172 
0173 g2   = cros;
0174 cros(:,1) = (cros(:,1)-ma)/sa;
0175 %cros0=cros;
0176 
0177 if 0, % slightly smoothing the crossing spectrum
0178   tmp=smooth(cros(1:end,1),cros(1:end,2),0.95,cros(1:end,1)); 
0179   ind=(tmp>0);
0180   cros(ind,2)=tmp; clear tmp ind
0181 end
0182 
0183 indz = (cros(:,2)==0);
0184 if any(indz)
0185   cros(indz,2) = inf; % this is done in order to avoid a warning message
0186 end
0187 cros(~indz,2) = -log(cros(~indz,2));
0188 
0189 
0190 % NB! the smooth function does not always extrapolate well outside the edges
0191 % causing poor estimate of g  
0192 % We may alleviate this problem by: forcing the extrapolation
0193 % to be linear outside the edges or choosing a lower value for csm1
0194 % or not to extrapolate at all in the first smoothing but instead
0195 % extrapolate in the second smoothing. (Possibly better since csm2<<csm1)
0196 % Therefore replacing the old call
0197 %scros=smooth(cros(10:ncr-10,1),cros(10:ncr-10,2),csm1,cros(1:ncr,1)); 
0198 % with
0199 inds = 1+Ne:ncr-Ne;% indices to points we are smoothing over
0200 inde = 1+Ne:ncr-Ne;% indices to points we are smoothing over and
0201              % possibly  extrapolating if length(inds)<length(inde)
0202 scros = smooth(cros(inds,1),cros(inds,2),csm1,cros(inde,1),def,cvar(inds)); 
0203 
0204 if plotflag>1
0205   %plot(cros(10:ncr-10,1),cros(10:ncr-10,2),'r')
0206   plot(cros(:,1),cros(:,2),'r'),
0207   hold on
0208   plot(cros(inde,1),scros,'b'),hold off
0209   title('First smoothing')
0210   ylabel('-log(cross)')
0211   xlabel('crossing level')
0212   legend('Not smoothed','Smoothed',0)
0213   %return
0214   pause(ptime)       
0215 end
0216 
0217 %  scros has to have a single minimum
0218 if 0,% old call
0219   [smin imin]=min(scros);%
0220 else % new call: checking if we have a single minimum 
0221  imin = findcross(diff(scros))+1; 
0222  smin = scros(imin); 
0223  if length(imin)~=1,
0224    disp(['Warning:  There are ' num2str(length(imin)) ' minima/' ...
0225        'maxima after the first smoothing'])  
0226    [smin ind] = min(smin);
0227    imin       = imin(ind);
0228  end
0229 end
0230 %imin=imin+30
0231 %smin = scros(imin)
0232 
0233 scros1 = sqrt(2*abs(scros-smin));
0234 %scros1(1:imin)=-scros1(1:imin);
0235 scros1(1:imin) = 2*scros1(imin)-scros1(1:imin);
0236 
0237 scros2 = smooth(cros(inde,1),scros1,csm2,uu,def,gvar(inde));
0238 
0239 g(:,2) = scros2';%*sa; %multiply with stdev 
0240 
0241 if nargout>2|plotflag>0,
0242   cros2 = cros; 
0243   [cmin icmin] = min(cros(:,2));
0244   cros2(:,2)   = sqrt(2*abs(cros(:,2)-cmin));
0245   cros2(1:icmin,2) = 2*cros2(icmin,2)-cros2(1:icmin,2); 
0246   g2(:,2)  = cros2(:,2);
0247 end
0248 
0249 if chkder~=0
0250    for ix = 1:5
0251     dy = diff(g(:,2));
0252     if any(dy<=0)
0253       warning('The empirical crossing spectrum is not sufficiently smoothed.')
0254       disp('        The estimated transfer function, g, is not ')
0255       disp('        a strictly increasing function.')
0256       dy(dy>0)=eps;
0257       gvar = -([dy;0]+[0;dy])/2+eps;
0258       g(:,2) = smooth(g(:,1),g(:,2),1,g(:,1),def,ix*gvar);
0259     else 
0260       break
0261     end
0262   end
0263 end
0264 % wrong!!!                
0265 %test=0.02*sqrt(sum((uu-scros2).^2));
0266 %                  5 
0267 %this is not sqrt(int (g(u)-u)^2 du)
0268 %                 -5
0269 % The correct is
0270 if 0,  %either
0271   test = sqrt((param(2)-param(1))/(param(3)-1)*sum((uu-scros2).^2));
0272 else % or
0273   %test=sqrt(simpson(uu,(uu-scros2).^2));
0274 % or
0275   test=sqrt(trapz(uu,(uu-scros2).^2));
0276 end 
0277 
0278 
0279 if plotflag>0,
0280   %% Plotchanges made by Joakim Elvander 970707.
0281   %% Plots will be initiated in the file funplot_1
0282   
0283   % %clf
0284   % plot(uu,scros2,'r')
0285   % hold on
0286   % plot(uu,uu,'g--')
0287   % pause
0288   
0289   
0290   if plotflag>1
0291     stairs(cros2(:,1),cros2(:,2),'b'),hold on
0292     plot(cros(inde,1),scros1,'r')
0293     plot(uu,scros2,'g--'),hold off
0294     title('Second smoothing')
0295     ylabel('Transfer function g(u)')
0296     xlabel('Crossing level u')
0297     legend('Not smoothed','Smoothed once','Smoothed twice',0)
0298     pause(ptime)
0299   end
0300   
0301   
0302   %stairs(cros2(:,1),cros2(:,2))  
0303   %axis([-5 5 -5 5])
0304   %axis('square')
0305   %hold off
0306  
0307   %% Temporarily
0308   trplot(g,g2,ma,sa)
0309   %  funplot_1(uu,scros2,cros2);
0310   %legend(['Smoothed (T='  num2str(test) ')'],'g(u)=u','Not smoothed',0)
0311   %ylabel('Transfer function g(u)')
0312   %xlabel('Crossing level u')
0313   
0314   if plotflag>1,pause(ptime),end
0315   %hold on, stairs(cros2(10:ncr-10,1),cros2(10:ncr-10,2),'y');hold off
0316 end
0317 
0318 
0319 
0320 
0321 
0322 
0323 
0324

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