LC2TR2 Estimate transformation, g, from observed crossing intensity, version2. Assumption: a Gaussian process, Y, is related to the non-Gaussian process, X, by Y = g(X). CALL: [g, test,g2] = lc2tr2(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 crossing intensity and the transformation g. 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 theoretical for a Gaussian model. 2 monitor development of 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 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) Ntr - Maximum length of empirical crossing intensity. The empirical crossing intensity is interpolated linearly before smoothing if the length 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) 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 GSM or set a larger variance for GVAR. 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 = lc2tr2(lc,0,Hm0/4,'plot','iter'); % Monitor the development g1 = lc2tr2(lc,0,Hm0/4,troptset('gvar', .5 )); % Equal weight on all points g2 = lc2tr2(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
Calculates discrete levels given the parameter matrix. | |
Calculates a smoothing spline. | |
Create or alter TRANSFORM OPTIONS structure. | |
Plots transformation, g, eg. estimated with dat2tr. | |
Inverse of the Normal distribution function | |
Cumulative trapezoidal numerical integration. | |
Difference and approximate derivative. | |
Display message and abort function. | |
1-D interpolation (table lookup) | |
Quick 1-D linear interpolation. | |
True if arrays are numerically equal. | |
Linearly spaced vector. | |
Wait for user response. | |
Trapezoidal numerical integration. | |
Display warning message; disable or enable warning messages. |
001 function [g, test, g2] = lc2tr(cross,ma,sa,varargin); 002 %LC2TR2 Estimate transformation, g, from observed crossing intensity, version2. 003 % 004 % Assumption: a Gaussian process, Y, is related to the 005 % non-Gaussian process, X, by Y = g(X). 006 % 007 % CALL: [g, test,g2] = lc2tr2(lc,ma,sa,options); 008 % 009 % g,g2 = smoothed and empirical estimate of the transformation g. 010 % test = test observator int (g(u)-u)^2 du where int limits is 011 % given by param. This is a measure of departure of the 012 % data from the Gaussian model. 013 % lc = a two column matrix with levels in the first column 014 % and number of upcrossings in the second. 015 % ma,sa = mean and standard deviation of the process 016 % 017 % options = structure with the fields: 018 % csm,gsm - defines the smoothing of the crossing intensity and the 019 % transformation g. Valid values must be 020 % 0<=csm,gsm<=1. (default csm = 0.9 gsm=0.05) 021 % Smaller values gives smoother functions. 022 % param - vector which defines the region of variation of the data X. 023 % (default [-5 5 513]). 024 % plotflag - 0 no plotting (Default) 025 % 1 plots empirical and smoothed g(u) and theoretical for a Gaussian model. 026 % 2 monitor development of estimation 027 % linextrap - 0 use a regular smoothing spline 028 % 1 use a smoothing spline with a constraint on the ends to 029 % ensure linear extrapolation outside the range of the data. 030 % (default) 031 % cvar - Variances for the crossing intensity. (default 1) 032 % gvar - Variances for the empirical transformation, g. (default 1) 033 % ne - Number of extremes (maxima & minima) to remove from the 034 % estimation of the transformation. This makes the 035 % estimation more robust against outliers. (default 7) 036 % Ntr - Maximum length of empirical crossing intensity. 037 % The empirical crossing intensity is interpolated 038 % linearly before smoothing if the length exceeds Ntr. 039 % A reasonable NTR will significantly speed up the 040 % estimation for long time series without loosing any 041 % accuracy. NTR should be chosen greater than 042 % PARAM(3). (default 1000) 043 % 044 % The empirical crossing intensity is usually very irregular. 045 % More than one local maximum of the empirical crossing intensity 046 % may cause poor fit of the transformation. In such case one 047 % should use a smaller value of GSM or set a larger variance for GVAR. 048 % If X(t) is likely to cross levels higher than 5 standard deviations 049 % then the vector param has to be modified. For example if X(t) is 050 % unlikely to cross a level of 7 standard deviations one can use 051 % param = [-7 7 513]. 052 % 053 % Example: 054 % Hm0 = 7; 055 % S = jonswap([],Hm0); g=ochitr([],[Hm0/4]); 056 % S.tr=g;S.tr(:,2)=g(:,2)*Hm0/4; 057 % xs = spec2sdat(S,2^13); 058 % lc = dat2lc(xs); 059 % g0 = lc2tr2(lc,0,Hm0/4,'plot','iter'); % Monitor the development 060 % g1 = lc2tr2(lc,0,Hm0/4,troptset('gvar', .5 )); % Equal weight on all points 061 % g2 = lc2tr2(lc,0,Hm0/4,'gvar', [3.5 .5 3.5]); % Less weight on the ends 062 % hold on, trplot(g1,g) % Check the fit 063 % trplot(g2) 064 % 065 % See also troptset, dat2tr, trplot, findcross, smooth 066 067 % NB! the transformated data will be N(0,1) 068 069 % Reference 070 % Rychlik , I., Johannesson, P., and Leadbetter, M.R. (1997) 071 % "Modelling and statistical analysis of ocean wavedata 072 % using a transformed Gaussian process", 073 % Marine structures, Design, Construction and Safety, 074 % Vol 10, pp 13--47 075 % 076 077 078 % Tested on: Matlab 5.3, 5.2, 5.1 079 % History: 080 % by pab 29.12.2000 081 % based on lc2tr, but the inversion is faster. 082 % by IR and PJ 083 084 085 opt = troptset('chkder','on','plotflag','off','csm',.9,'gsm',.05,.... 086 'param',[-5 5 513],'delay',2,'linextrap','on','ntr',1000,'ne',7,'cvar',1,'gvar',1); 087 % If just 'defaults' passed in, return the default options in g 088 if nargin==1 & nargout <= 1 & isequal(cross,'defaults') 089 g = opt; 090 return 091 end 092 error(nargchk(3,inf,nargin)) 093 if nargin>=4 , opt = troptset(opt,varargin{:}); end 094 csm2 = opt.gsm; 095 param = opt.param; 096 ptime = opt.delay; 097 Ne = opt.ne; 098 switch opt.chkder; 099 case 'off', chkder = 0; 100 case 'on', chkder = 1; 101 otherwise, chkder = opt.chkder; 102 end 103 switch opt.linextrap; 104 case 'off', def = 0; 105 case 'on', def = 1; 106 otherwise, def = opt.linextrap; 107 end 108 switch opt.plotflag 109 case {'none','off'}, plotflag = 0; 110 case 'final', plotflag = 1; 111 case 'iter', plotflag = 2; 112 otherwise, plotflag = opt.plotflag; 113 end 114 ncr = length(cross); 115 if ncr>opt.ntr & opt.ntr>0, 116 x0 = linspace(cross(1+Ne,1),cross(end-Ne,1),opt.ntr)'; 117 cros = [ x0,interp1q(cross(:,1),cross(:,2),x0)]; 118 Ne = 0; 119 Ner = opt.ne; 120 ncr = opt.ntr; 121 else 122 Ner = 0; 123 cros=cross; 124 end 125 126 ng = length(opt.gvar); 127 if ng==1 128 gvar = opt.gvar(ones(ncr,1)); 129 else 130 gvar = interp1(linspace(0,1,ng)',opt.gvar(:),linspace(0,1,ncr)','*linear'); 131 end 132 ng = length(opt.cvar); 133 if ng==1 134 cvar = opt.cvar(ones(ncr,1)); 135 else 136 cvar = interp1(linspace(0,1,ng)',opt.cvar(:),linspace(0,1,ncr)','*linear'); 137 end 138 139 g = zeros(param(3),2); 140 141 uu = levels(param); 142 143 g(:,1) = sa*uu' + ma; 144 145 g2 = cros; 146 147 if Ner>0, % Compute correction factors 148 cor1 = trapz(cross(1:Ner+1,1),cross(1:Ner+1,2)); 149 cor2 = trapz(cross(end-Ner-1:end,1),cross(end-Ner-1:end,2)); 150 else 151 cor1 = 0; 152 cor2 = 0; 153 end 154 cros(:,2) = cumtrapz(cros(:,1),cros(:,2))+cor1; 155 cros(:,2) = (cros(:,2)+.5)/(cros(end,2) + cor2 +1); 156 cros(:,1) = (cros(:,1)-ma)/sa; 157 158 % find the mode 159 [tmp,imin]= min(abs(cros(:,2)-.15)); 160 [tmp,imax]= min(abs(cros(:,2)-.85)); 161 inde = imin:imax; 162 tmp = smooth(cros(inde,1),g2(inde,2),opt.csm,cros(inde,1),def,cvar(inde)); 163 164 [tmp imax] = max(tmp); 165 u0 = cros(inde(imax),1); 166 %u0 = interp1q(cros(:,2),cros(:,1),.5) 167 168 169 cros(:,2) = wnorminv(cros(:,2),-u0,1); 170 171 g2(:,2) = cros(:,2); 172 % NB! the smooth function does not always extrapolate well outside the edges 173 % causing poor estimate of g 174 % We may alleviate this problem by: forcing the extrapolation 175 % to be linear outside the edges or choosing a lower value for csm2. 176 177 inds = 1+Ne:ncr-Ne;% indices to points we are smoothing over 178 scros2 = smooth(cros(inds,1),cros(inds,2),csm2,uu,def,gvar(inds)); 179 180 g(:,2) = scros2';%*sa; %multiply with stdev 181 182 if chkder~=0 183 for ix = 1:5 184 dy = diff(g(:,2)); 185 if any(dy<=0) 186 warning('The empirical crossing spectrum is not sufficiently smoothed.') 187 disp(' The estimated transfer function, g, is not ') 188 disp(' a strictly increasing function.') 189 dy(dy>0)=eps; 190 gvar = -([dy;0]+[0;dy])/2+eps; 191 g(:,2) = smooth(g(:,1),g(:,2),1,g(:,1),def,ix*gvar); 192 else 193 break 194 end 195 end 196 end 197 if 0, %either 198 test = sqrt((param(2)-param(1))/(param(3)-1)*sum((uu-scros2).^2)); 199 else % or 200 %test=sqrt(simpson(uu,(uu-scros2).^2)); 201 % or 202 test=sqrt(trapz(uu,(uu-scros2).^2)); 203 end 204 205 206 if plotflag>0, 207 trplot(g ,g2,ma,sa) 208 %legend(['Smoothed (T=' num2str(test) ')'],'g(u)=u','Not smoothed',0) 209 %ylabel('Transfer function g(u)') 210 %xlabel('Crossing level u') 211 212 if plotflag>1,pause(ptime),end 213 end 214 215 216 217 218 219 220 221
Comments or corrections to the WAFO group