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)

## CROSS-REFERENCE INFORMATION

This function calls:
 cdf2tr Estimate transformation, g, from observed CDF. dat2tp Extracts turning points from data, empdistr Computes and plots the empirical CDF hermitetr Calculate transformation, g, proposed by Winterstein lc2tr Estimate transformation, g, from observed crossing intensity. levels Calculates discrete levels given the parameter matrix. mm2lc Extracts level-crossing spectrum from min2Max cycles. ochitr Calculates transformation, g, proposed by Ochi et al. smooth Calculates a smoothing spline. tp2mm Calculates min2Max and Max2min cycles from a sequence of turning points troptset Create or alter TRANSFORM OPTIONS structure. wkurtosis Computes sample kurtosis wskewness Computes sample skewness error Display message and abort function. isequal True if arrays are numerically equal. lower Convert string to lowercase. mean Average or mean value. pause Wait for user response. std Standard deviation. strcmpi Compare strings ignoring case. trapz Trapezoidal numerical integration.
This function is called by:
 Chapter2 % CHAPTER2 Modelling random loads and stochastic waves Chapter3 % CHAPTER3 Demonstrates distributions of wave characteristics dat2spec Estimate one-sided spectral density from data. dat2spec2 Estimate one-sided spectral density, version 2. mctrtest Test if a stochastic process is Gaussian. reconstruct reconstruct the spurious/missing points of timeseries

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