CDF2TR Estimate transformation, g, from observed CDF. Assumption: a Gaussian process, Y, is related to the non-Gaussian process, X, by Y = g(X). CALL [g,test,g2] = cdf2tr(F,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 OPTIONS.PARAM. This is a measure of departure of the data from the Gaussian model. F = empirical CDF of X(t), a 2 column matrix. ma,sa = mean and standard deviation of the process X(t). options = options structure defining how the smoothing is done. (See troptset for default values) The empirical CDF is usually very irregular. More than one local maximum of the empirical CDF 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); Fx = empdistr(xs(:,2)); g0 = cdf2tr(Fx,0,Hm0/4,troptset('plot',1)); % Plot final estimate g1 = lc2tr(Fx,0,Hm0/4,troptset('gvar', .5 )); % More weight on all points g2 = lc2tr(Fx,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, lc2tr
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 | |
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. |
Estimate transformation, g, from data. |
001 function [g, test, g2] = cdf2tr(Fx1,ma ,sa,varargin) 002 %CDF2TR Estimate transformation, g, from observed CDF. 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] = cdf2tr(F,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 OPTIONS.PARAM. This is a measure of departure of the 012 % data from the Gaussian model. 013 % F = empirical CDF of X(t), a 2 column matrix. 014 % ma,sa = mean and standard deviation of the process X(t). 015 % options = options structure defining how the smoothing is done. 016 % (See troptset for default values) 017 % 018 % The empirical CDF is usually very irregular. 019 % More than one local maximum of the empirical CDF 020 % may cause poor fit of the transformation. In such case one 021 % should use a smaller value of GSM or set a larger variance for GVAR. 022 % If X(t) is likely to cross levels higher than 5 standard deviations 023 % then the vector param has to be modified. For example if X(t) is 024 % unlikely to cross a level of 7 standard deviations one can use 025 % param = [-7 7 513]. 026 % 027 % Example 028 % Hm0 = 7; 029 % S = jonswap([],Hm0); g=ochitr([],[Hm0/4]); 030 % S.tr = g; S.tr(:,2)=g(:,2)*Hm0/4; 031 % xs = spec2sdat(S,2^13); 032 % Fx = empdistr(xs(:,2)); 033 % g0 = cdf2tr(Fx,0,Hm0/4,troptset('plot',1)); % Plot final estimate 034 % g1 = lc2tr(Fx,0,Hm0/4,troptset('gvar', .5 )); % More weight on all points 035 % g2 = lc2tr(Fx,0,Hm0/4,'gvar', [3.5 .5 3.5]); % Less weight on the ends 036 % hold on, trplot(g1,g) % Check the fit 037 % trplot(g2) 038 % 039 % See also troptset, lc2tr 040 041 % History 042 % revised Feb2004 043 % Revised pab Dec2003 044 % fixed a bug F -> Fx1 045 % by pab 29.12.2000 046 % - default param is now [-5 5 513] -> better to have the discretization 047 % represented with exact numbers, especially when calculating 048 % derivatives of the transformation numerically. 049 050 opt = troptset('chkder','on','plotflag','off','gsm',.05,.... 051 'param',[-5 5 513],'delay',2,'linextrap','on','ntr',1000,'ne',7,'gvar',1); 052 % If just 'defaults' passed in, return the default options in g 053 if nargin==1 & nargout <= 1 & isequal(Fx1,'defaults') 054 g = opt; 055 return 056 end 057 error(nargchk(3,inf,nargin)) 058 if nargin>=4, opt=troptset(opt,varargin{:}); end 059 switch opt.chkder; 060 case 'off', chkder = 0; 061 case 'on', chkder = 1; 062 otherwise, chkder = opt.chkder; 063 end 064 switch opt.linextrap; 065 case 'off', def = 0; 066 case 'on', def = 1; 067 otherwise, def = opt.linextrap; 068 end 069 070 switch opt.plotflag 071 case {'none','off'}, plotflag = 0; 072 case 'final', plotflag = 1; 073 case 'iter', plotflag = 2; 074 otherwise, plotflag = opt.plotflag; 075 end 076 077 078 Ne = opt.ne; 079 if length(Fx1)>opt.ntr & opt.ntr>0 080 x0 = linspace(Fx1(1+Ne,1),Fx1(end-Ne,1),opt.ntr)'; 081 Fx = [ x0,interp1q(Fx1(:,1),Fx1(:,2),x0)]; 082 Ne=0; 083 else 084 Fx = Fx1; 085 end 086 uu = levels(opt.param)'; 087 g = [sa*uu+ma zeros(opt.param(3),1)]; 088 ncr = length(Fx); 089 090 091 ng = length(opt.gvar); 092 if ng==1 093 gvar = opt.gvar(ones(ncr,1)); 094 else 095 gvar = interp1(linspace(0,1,ng)',opt.gvar(:),linspace(0,1,ncr)','*linear'); 096 end 097 098 099 ind = find(diff(Fx(:,1))>0); % remove equal points 100 ind1 = ind(Ne+1:end-Ne); 101 tmp = wnorminv(Fx(ind,2)); 102 103 g(:,2) = smooth(Fx(ind1,1),tmp(Ne+1:end-Ne),opt.gsm,g(:,1),def,gvar); 104 105 if chkder~=0 106 for ix = 1:5 107 dy = diff(g(:,2)); 108 if any(dy<=0) 109 warning('The empirical distribution is not sufficiently smoothed.') 110 disp(' The estimated transfer function, g, is not ') 111 disp(' a strictly increasing function.') 112 dy(dy>0)=eps; 113 gvar = -([dy;0]+[0;dy])/2+eps; 114 g(:,2) = smooth(g(:,1),g(:,2),1,g(:,1),def,ix*gvar); 115 else 116 break 117 end 118 end 119 end 120 if nargout>1 121 test = sqrt(trapz(uu,(uu-g(:,2)).^2)); 122 end 123 if nargout>2 124 g2 = [Fx(ind,1) tmp]; 125 end 126 if plotflag>0 127 trplot(g,g2,ma,sa) 128 if plotflag>1 129 pause(opt.delay) 130 end 131 end 132 133 134 135
Comments or corrections to the WAFO group