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)

## CROSS-REFERENCE INFORMATION

This function calls:
 findcross Finds indices to level v up and downcrossings of a vector levels Calculates discrete levels given the parameter matrix. smooth Calculates a smoothing spline. troptset Create or alter TRANSFORM OPTIONS structure. trplot Plots transformation, g, eg. estimated with dat2tr. clear Clear variables and functions from memory. diff Difference and approximate derivative. error Display message and abort function. hold Hold current graph. interp1 1-D interpolation (table lookup) interp1q Quick 1-D linear interpolation. isequal True if arrays are numerically equal. legend Display legend. linspace Linearly spaced vector. num2str Convert number to string. (Fast version) pause Wait for user response. plot Linear plot. stairs Stairstep plot. title Graph title. trapz Trapezoidal numerical integration. warning Display warning message; disable or enable warning messages. xlabel X-axis label. ylabel Y-axis label.
This function is called by:
 dat2tr Estimate transformation, g, from data. lc2sdat Simulates process with given irregularity factor and crossing spectrum

## 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 %
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