Home > wafo > trgauss > dat2tr.m

dat2tr

PURPOSE ^

Estimate transformation, g, from data.

SYNOPSIS ^

[g, test, cmax, irr, g2]= dat2tr(x,def,varargin);

DESCRIPTION ^

 DAT2TR Estimate transformation, g, from data.
 
  CALL:  [g test cmax irr g2]  = dat2tr(x,def,options);
 
    g,g2   = the smoothed and empirical transformation, respectively. 
             A two column matrix if multip=0.  
             If multip=1 it ís a 2*(m-1) column matrix where the
             first and second column is the transform 
             for values in column 2 and third and fourth column is the
             transform for values in column 3 ......
 
    test   = 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.
            
    cmax   = maximum crossing intensity of x
    irr    = irregularity factor of x which is approximately Tz/Tmaxima   
    x      = m column data matrix with sampled times in the first column
             and values the next columns.            
 
    def    = 'nonlinear' : transform based on smoothed crossing intensity (default)
             'mnonlinear': transform based on smoothed marginal distribution
             'hermite'   : transform based on cubic Hermite polynomial
             'ochitr'    : transform based on exponential function
             'linear'    : identity.
 
    options = options structure with the following 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 see lc2tr). 
  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)
      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 or CDF.
             The empirical crossing intensity or CDF is interpolated
             linearly  before smoothing if their lengths 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)
    multip - 0 the data in columns belong to the same seastate (default).
             1 the data in columns are from separate seastates.
 
   DAT2TR estimates the transformation in a transformed Gaussian model.  
   Assumption: a Gaussian process, Y, is related to the
   non-Gaussian process, X, by Y = g(X). 
  
   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 CSM. In order to check the effect 
   of smoothing it is recomended to also plot g and g2 in the same plot or
   plot the smoothed g against an interpolated version of g (when CSM=GSM=1).
     If  x  is likely to cross levels higher than 5 standard deviations
   then the vector param has to be modified.  For example if x 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);
  g0 = dat2tr(xs,[],'plot','iter');             % Monitor the development
  g1 = dat2tr(xs,'mnon','gvar', .5 );           % More weight on all points
  g2 = dat2tr(xs,'nonl','gvar', [3.5 .5 3.5]);  % Less weight on the ends
  hold on, trplot(g1,g)                                   % Check the fit
  trplot(g2)
 
  See also  troptset, lc2tr, cdf2tr, trplot

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [g, test, cmax, irr, g2]= dat2tr(x,def,varargin);
002 %DAT2TR Estimate transformation, g, from data.
003 %
004 % CALL:  [g test cmax irr g2]  = dat2tr(x,def,options);
005 %
006 %   g,g2   = the smoothed and empirical transformation, respectively. 
007 %            A two column matrix if multip=0.  
008 %            If multip=1 it ís a 2*(m-1) column matrix where the
009 %            first and second column is the transform 
010 %            for values in column 2 and third and fourth column is the
011 %            transform for values in column 3 ......
012 %
013 %   test   = int (g(u)-u)^2 du  where int. limits is given by param. This
014 %            is a measure of departure of the data from the Gaussian model.
015 %           
016 %   cmax   = maximum crossing intensity of x
017 %   irr    = irregularity factor of x which is approximately Tz/Tmaxima   
018 %   x      = m column data matrix with sampled times in the first column
019 %            and values the next columns.            
020 %
021 %   def    = 'nonlinear' : transform based on smoothed crossing intensity (default)
022 %            'mnonlinear': transform based on smoothed marginal distribution
023 %            'hermite'   : transform based on cubic Hermite polynomial
024 %            'ochitr'    : transform based on exponential function
025 %            'linear'    : identity.
026 %
027 %   options = options structure with the following fields:
028 %  csm,gsm - defines the smoothing of the logarithm of crossing intensity 
029 %            and the transformation g, respectively. Valid values must 
030 %            be 0<=csm,gsm<=1. (default csm=0.9, gsm=0.05)
031 %            Smaller values gives smoother functions.
032 %    param - vector which defines the region of variation of the data x.
033 %           (default see lc2tr). 
034 % plotflag - 0 no plotting (Default)
035 %            1 plots empirical and smoothed g(u) and the theoretical for
036 %              a Gaussian model. 
037 %            2 monitor the development of the estimation
038 %linextrap - 0 use a regular smoothing spline 
039 %            1 use a smoothing spline with a constraint on the ends to 
040 %              ensure linear extrapolation outside the range of the data.
041 %              (default)
042 %     gvar - Variances for the empirical transformation, g. (default  1) 
043 %       ne - Number of extremes (maxima & minima) to remove from the
044 %            estimation of the transformation. This makes the
045 %            estimation more robust against outliers. (default 7)
046 %      ntr - Maximum length of empirical crossing intensity or CDF.
047 %            The empirical crossing intensity or CDF is interpolated
048 %            linearly  before smoothing if their lengths exceeds Ntr.
049 %            A reasonable NTR will significantly speed up the
050 %            estimation for long time series without loosing any
051 %            accuracy. NTR should be chosen greater than
052 %            PARAM(3). (default 1000)
053 %   multip - 0 the data in columns belong to the same seastate (default).
054 %            1 the data in columns are from separate seastates.
055 %
056 %  DAT2TR estimates the transformation in a transformed Gaussian model.  
057 %  Assumption: a Gaussian process, Y, is related to the
058 %  non-Gaussian process, X, by Y = g(X). 
059 % 
060 %  The empirical crossing intensity is usually very irregular.
061 %  More than one local maximum of the empirical crossing intensity
062 %  may cause poor fit of the transformation. In such case one
063 %  should use a smaller value of CSM. In order to check the effect 
064 %  of smoothing it is recomended to also plot g and g2 in the same plot or
065 %  plot the smoothed g against an interpolated version of g (when CSM=GSM=1).
066 %    If  x  is likely to cross levels higher than 5 standard deviations
067 %  then the vector param has to be modified.  For example if x is 
068 %  unlikely to cross a level of 7 standard deviations one can use 
069 %  PARAM=[-7 7 513].
070 %
071 % Example:
072 % Hm0 = 7;
073 % S = jonswap([],Hm0); g=ochitr([],[Hm0/4]); 
074 % S.tr=g;S.tr(:,2)=g(:,2)*Hm0/4;
075 % xs = spec2sdat(S,2^13);
076 % g0 = dat2tr(xs,[],'plot','iter');             % Monitor the development
077 % g1 = dat2tr(xs,'mnon','gvar', .5 );           % More weight on all points
078 % g2 = dat2tr(xs,'nonl','gvar', [3.5 .5 3.5]);  % Less weight on the ends
079 % hold on, trplot(g1,g)                                   % Check the fit
080 % trplot(g2)
081 %
082 % See also  troptset, lc2tr, cdf2tr, trplot
083 
084 % References:
085 % Rychlik, I. , Johannesson, P and Leadbetter, M. R. (1997)
086 % "Modelling and statistical analysis of ocean wavedata using 
087 %  transformed Gaussian process."
088 % Marine structures, Design, Construction and Safety, Vol. 10, No. 1, pp 13--47
089 %
090 % 
091 % Brodtkorb, P, Myrhaug, D, and Rue, H (1999)
092 % "Joint distribution of wave height and crest velocity from
093 % reconstructed data"
094 % in Proceedings of 9th ISOPE Conference, Vol III, pp 66-73
095 
096 
097 
098 %Tested on: Matlab 5.3, 5.2, 5.1
099 %History:
100 % revised pab Dec2004
101 %  -Fixed bug: string comparison for def at fault.  
102 % revised pab Nov2004
103 %  -Fixed bug: linextrap was not accounted for  
104 % revised pab july 2004
105 % revised pab 3 april 2004
106 % -fixed a bug in hermite estimation: excess changed to kurtosis  
107 % revised pab 29.12.2000
108 % - added example, hermite and ochi options
109 % - replaced optional arguments with a options struct
110 % - default param is now [-5 5 513] -> better to have the discretization
111 %  represented with exact numbers, especially when calculating
112 %  derivatives of the transformation numerically.
113 % revised pab 19.12.2000
114 %  - updated call empdistr(X,-inf,[],monitor) to  empdistr(X,[],monitor)
115 %    due to new calling syntax for empdistr
116 % modifed pab 24.09.2000
117 %  - changed call from norminv to wnorminv
118 %  - also removed the 7 lowest and 7 highest points from
119 %    the estimation using def='mnonlinear' 
120 %    (This is similar to what lc2tr does. lc2tr removes
121 %     the 9 highest and 9 lowest TP from the estimation)
122 % modified pab 09.06.2000
123 %  - made all the *empirical options secret.
124 %  - Added 'mnonlinear' and 'mempirical' 
125 %  - Fixed the problem of multip==1 and def=='empirical' by interpolating 
126 %    with spline to ensure that the length of g is fixed
127 %  - Replaced the test statistic for def=='empirical' with the one
128 %    obtained when csm1=csm2=1. (Previously only the smoothed test
129 %    statistic where returned)
130 % modified pab 12.10.1999
131 %  fixed a bug
132 %  added secret output of empirical estimate g2
133 % modified by svi  29.09.1999
134 % changed input def by adding new options.
135 % revised by pab 11.08.99
136 %   changed name from dat2tran to dat2tr
137 % modified by Per A. Brodtkorb 12.05.1999,15.08.98
138 %   added  secret option: to accept multiple data, to monitor the steps 
139 %   of estimation of the transformation 
140 %   also removed some code and replaced it with a call to lc2tr (cross2tr) 
141 %   making the maintainance easier
142 %
143 
144 %opt = troptset('plotflag','off','csm',.95,'gsm',.05,....
145 %    'param',[-5 5 513],'delay',2,'linextrap','on','ne',7,...
146 %    'cvar',1,'gvar',1,'multip',0);
147 
148 
149 opt = troptset('chkder','on','plotflag','off','csm',.95,'gsm',.05,....
150     'param',[-5 5 513],'delay',2,'ntr',1000,'linextrap','on','ne',7,'cvar',1,'gvar',1,'multip',0,'crossdef','uM');
151 % If just 'defaults' passed in, return the default options in g
152 if nargin==1 & nargout <= 1 & isequal(x,'defaults')
153   g = opt; 
154   return
155 end
156 error(nargchk(1,inf,nargin)) 
157 if nargin<2|isempty(def),     def    = 'nonlinear'; end
158 if nargin>=3,  opt   = troptset(opt,varargin{:});   end
159 multip = opt.multip;
160 Ne = opt.ne;
161 switch opt.plotflag
162   case {'none','off'},   plotflag = 0;
163   case 'final', plotflag = 1;
164  case 'iter',  plotflag = 2;
165   
166   otherwise,    plotflag = opt.plotflag;
167 end
168 
169 % Crossing definition
170 switch opt.crossdef
171   case {'u'}, cdef = 1; % only upcrossings.
172   case 'uM',  cdef = 2; % upcrossings and Maxima (default).
173   case 'umM', cdef = 3; % upcrossings, minima, and Maxima.
174   case 'um',  cdef = 4; % upcrossings and minima.  
175   otherwise,  cdef = opt.crossdef;
176 end
177 switch opt.linextrap
178  case {'on'},  linextrap = 1;
179  case {'off'}, linextrap = 0;
180  otherwise,    linextrap = opt.linextrap;
181 end
182 
183 xx = x;
184 [n,m] = size(xx);
185 ma = mean(xx(:,2:m));
186 sa = std(xx(:,2:m)); 
187 m2 = m;
188 
189 if m>2 & multip==0,% data in columns belongs to the same seastate
190   ma = mean(ma);
191   sa = sqrt(mean(sa.^2)); % pooled standard deviation
192   m2 = 2;
193 end
194 
195   
196 g  = zeros(opt.param(3),2*(m2-1));
197 uu = levels(opt.param)';
198 g(:,2:2:end) = uu(:,ones(1,m2-1));
199 g(:,1:2:end) = sa(ones(opt.param(3),1),:).*g(:,2:2:end) + ...
200       ma(ones(opt.param(3),1),:);
201 g2 = g;
202 
203 test = zeros(m2-1,1);
204 
205 if strcmpi(lower(def(1:3)),'lin') & nargout<=2, return,end
206 
207 irr   = test;
208 cmax  = irr;
209 
210 
211 if multip==1,
212   for ix=1:(m-1),
213     if (lower(def(1))=='n'|lower(def(1))=='e' | nargout>2),
214       tp       = dat2tp(xx(:,[1 ix+1]));% Turning points 
215       mM       = tp2mm(tp);             % min2Max cycles
216       cross1   = mm2lc(mM,cdef,0);      % Want upcrossings and maxima
217       cmax(ix) = max(cross1(:,2));
218       irr(ix)  = length(mM)/cmax(ix);   % approximately Tz/Tmaxima
219     end
220     if (lower(def(1))=='m')
221       Fx  = empdistr(xx(:,[ ix+1]),[],plotflag==2);
222       if plotflag==2
223     pause(opt.delay)
224       end
225     end
226     switch lower(def(1:3))
227       case {'her'},
228     ga1 = wskewness(xx(:,ix+1));
229     ga2 = wkurtosis(xx(:,ix+1))-3;
230     up  = min(4*(4*ga1/3).^2,12);
231     lo  = ga1^2*3/2;
232     kurt = min(up,max(ga2,lo))+3;
233     phat = [sa(ix), ga1, kurt,ma(ix) ];
234     [g(:,2*ix-1:2*ix), test(ix)] = hermitetr(g(:,2*ix-1) ,phat);
235     g2=[];
236       case {'och'},
237     phat = [ sa(ix) wskewness(xx(:,ix+1)) ma(ix)];
238     [g(:,2*ix-1:2*ix), test(ix)] = ochitr(g(:,2*ix-1) ,phat);
239     g2=[];
240       case {'non'}, % nonlinear      
241     [g(:,2*ix-1:2*ix), test(ix),tmp]=lc2tr(cross1,ma(ix),sa(ix),opt);
242     if nargout>4
243       g2(:,2*ix-1) = g(:,2*ix-1);
244       g2(:,2*ix)   = smooth(tmp(:,1),tmp(:,2),1,g(:,2*ix-1),linextrap);
245     end
246       case {'emp'}, % empirical 
247     [g2(:,2*ix-1:2*ix), test(ix),tmp]=lc2tr(cross1,ma(ix),sa(ix),opt);
248     g(:,2*ix-1) = g2(:,2*ix-1);
249     g(:,2*ix)   = smooth(tmp(:,1),tmp(:,2),1,g2(:,2*ix-1),1);
250     test(ix)    = sqrt(trapz(uu,(uu-g(:,2*ix)).^2));
251       case {'mno'} % mnonlinear 
252     [g(:,2*ix-1:2*ix), test(ix),tmp]= cdf2tr(Fx,ma(ix),sa(ix),opt);
253     if nargout>4
254       g2(:,2*ix-1) = g(:,2*ix-1);
255       g2(:,2*ix)   = smooth(tmp(:,1),tmp(:,2),1,g(:,2*ix-1),1);
256     end
257       case {'mem'}, % mempirical
258     [g2(:,2*ix-1:2*ix), test(ix),tmp]=cdf2tr(Fx,ma(ix),sa(ix),opt);
259     g(:,2*ix-1) = g2(:,2*ix-1);
260     g(:,2*ix)   = smooth(tmp(:,1),tmp(:,2),1,g2(:,2*ix-1),linextrap);
261     test(ix)    = sqrt(trapz(uu,(uu-g(:,2*ix)).^2));
262     %g(:,2*ix)  = smooth(Fx(ind,1),tmp,1,g(:,2*ix-1),linextrap);
263     %test(ix)   = sqrt(trapz(uu,(uu-g(:,2*ix)).^2));
264     end
265   end
266 else % multip==0
267   if (lower(def(1))=='n'|lower(def(1))=='e' | nargout>2),
268     tp=[];mM=[];
269     for ix=1:(m-1),
270       tmp=dat2tp(xx(:,[1 ix+1]));
271       tp=[tp; tmp];
272       mM=[mM;tp2mm(tmp)];
273     end
274     
275     cross1 = mm2lc(mM,cdef,0); %want upcrossings and maxima  
276     cmax   = max(cross1(:,2));
277     irr    = length(mM)/cmax;% approximately Tz/Tmaxima
278   end
279   if (lower(def(1))=='m')
280     Fx  = empdistr(xx(n+1:end),[],plotflag==2);
281     if plotflag==2
282       pause(opt.delay)
283     end
284   end
285   switch lower(def(1:3))
286     case {'her'},
287       ga1 = wskewness(xx(n+1:end));
288       ga2 = wkurtosis(xx(n+1:end))-3;
289       up  = min(4*(4*ga1/3).^2,13);
290       lo  = ga1^2*3/2;
291       kurt = min(up,max(ga2,lo))+3;
292       phat = [sa ga1,kurt,ma ];
293       [g, test] = hermitetr(g(:,1) ,phat);
294       g2=[];
295     case {'och'}
296       phat = [sa wskewness(xx(n+1:end)) ma];
297       [g test] = ochitr(g(:,1) ,phat);
298       g2=[];
299     case {'non'},
300       [g, test, g2] = lc2tr(cross1,ma,sa,opt);%csm1,csm2,param,[plotflag monitor]);
301     case {'emp'},  
302       [g2, test, g] = lc2tr(cross1,ma,sa,opt);%csm1,csm2,param,[plotflag monitor]);
303       test = sqrt(trapz(uu,(uu-smooth((g(:,1)-ma)/sa,g(:,2),1,uu,1)).^2));
304     case {'mno'},
305       [g, test, g2] = cdf2tr(Fx,ma,sa,opt);
306     case  {'mem'},
307       [g2, test, g] = cdf2tr(Fx,ma,sa,opt);
308       test = sqrt(trapz(uu,(uu-smooth((g(:,1)-ma)/sa,g(:,2),1,uu,1)).^2));
309   end 
310 end
311

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