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
Estimate transformation, g, from observed CDF. | |
Extracts turning points from data, | |
Computes and plots the empirical CDF | |
Calculate transformation, g, proposed by Winterstein | |
Estimate transformation, g, from observed crossing intensity. | |
Calculates discrete levels given the parameter matrix. | |
Extracts level-crossing spectrum from min2Max cycles. | |
Calculates transformation, g, proposed by Ochi et al. | |
Calculates a smoothing spline. | |
Calculates min2Max and Max2min cycles from a sequence of turning points | |
Create or alter TRANSFORM OPTIONS structure. | |
Computes sample kurtosis | |
Computes sample skewness | |
Display message and abort function. | |
True if arrays are numerically equal. | |
Convert string to lowercase. | |
Average or mean value. | |
Wait for user response. | |
Standard deviation. | |
Compare strings ignoring case. | |
Trapezoidal numerical integration. |
% CHAPTER2 Modelling random loads and stochastic waves | |
% CHAPTER3 Demonstrates distributions of wave characteristics | |
Estimate one-sided spectral density from data. | |
Estimate one-sided spectral density, version 2. | |
Test if a stochastic process is Gaussian. | |
reconstruct the spurious/missing points of timeseries |
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
Comments or corrections to the WAFO group