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) See also troptset, dat2tr, trplot, findcross, smooth
Finds indices to level v up and downcrossings of a vector | |
Calculates discrete levels given the parameter matrix. | |
Calculates a smoothing spline. | |
Create or alter TRANSFORM OPTIONS structure. | |
Plots transformation, g, eg. estimated with dat2tr. | |
Clear variables and functions from memory. | |
Difference and approximate derivative. | |
Display message and abort function. | |
Hold current graph. | |
1-D interpolation (table lookup) | |
Quick 1-D linear interpolation. | |
True if arrays are numerically equal. | |
Display legend. | |
Linearly spaced vector. | |
Convert number to string. (Fast version) | |
Wait for user response. | |
Linear plot. | |
Stairstep plot. | |
Graph title. | |
Trapezoidal numerical integration. | |
Display warning message; disable or enable warning messages. | |
X-axis label. | |
Y-axis label. |
Estimate transformation, g, from data. | |
Simulates process with given irregularity factor and crossing spectrum |
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 % 0058 % See also troptset, dat2tr, trplot, findcross, smooth 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
Comments or corrections to the WAFO group