OCHITR Calculates transformation, g, proposed by Ochi et al. Assumption: a Gaussian process, Y, is related to the non-Gaussian process, X, by Y = g(X). CALL: [g, test]= ochitr(x,data); g = [x g(x)] a two column matrix with the transformation g(x). test = int (g(x)-x)^2 dy where int. limits is given by Y. This is a measure of departure of the data from the Gaussian model. x = vector with x-values. (default linspace(-5*sigma,5*sigma,513)+mean) data = the data vector [sigma skew mean] is the standard deviation, skewness and mean of the process, respectively. skew=0 for a Gaussian process. (abs(skew) <=2.82)(default [1 0.16 0]) This is a transformation model where the transformation is chosen to be a monotonic exponential function, calibrated such that the first 3 moments of the transformed model G(y)=g^-1(y) match the moments of the true process. However, the skewness is limited by ABS(SKEW)<2.82. Information about the moments of the process can be obtained by site specific data, laboratory measurements or by resort to theoretical models (see spec2skew). According to Ochi it is appropriate for a process with very strong non-linear characteristics. g(x) = ((1-exp(-gamma*(x-mean)/sigma))/gamma-mean2)/sigma2 where gamma = 1.28*a for x>=mean 3*a otherwise mean, sigma = standard deviation and mean, respectively, of the process. mean2, sigma2 = normalizing parameters in the transformed world, i.e., to make the gaussian process in the transformed world is N(0,1). The unknown parameters a, mean2 and sigma2 are found by solving the following non-linear equations: a*(sigma2^2+mean2^2)+mean2 = 0 sigma2^2-2*a^2*sigma2^4 = 1 2*a*sigma2^4*(3-8*a^2*sigma2^2) = skew NOTE: - by specifying NaN's in the data vector default values will be used. - if length(data) is shorter than the parameters needed then the default values are used for the parameters not specified. - The gaussian process in the transformed world is N(0,1) - g does not have continous derivatives of 2'nd order or higher. Example: Simulate a Transformed Gaussian process: Hm0=7;Tp=11; S = jonswap([],[Hm0 Tp]); [sk ku ma]=spec2skew(S); g = ochitr([],[Hm0/4,sk,ma]); g2=[g(:,1), g(:,2)*Hm0/4]; ys = spec2sdat(S,15000); % Simulated in the Gaussian world xs = gaus2dat(ys,g2); % Transformed to the real world See also spec2skew, hermitetr, lc2tr, dat2tr
Calculates the transformation g proposed by Ochi et al. | |
Display message and abort function. | |
Multidimensional unconstrained nonlinear minimization (Nelder-Mead). | |
Scalar nonlinear zero finding. | |
Linearly spaced vector. | |
Convert string matrix to numeric array. | |
MATLAB version number. |
Estimate transformation, g, from data. |
001 function [g,t0]=ochitr(y,data) 002 %OCHITR Calculates transformation, g, proposed by Ochi et al. 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]= ochitr(x,data); 008 % 009 % g = [x g(x)] a two column matrix with the transformation g(x). 010 % test = int (g(x)-x)^2 dy where int. limits is given by Y. This 011 % is a measure of departure of the data from the Gaussian model. 012 % x = vector with x-values. 013 % (default linspace(-5*sigma,5*sigma,513)+mean) 014 % data = the data vector [sigma skew mean] 015 % is the standard deviation, skewness and mean 016 % of the process, respectively. skew=0 for a Gaussian process. 017 % (abs(skew) <=2.82)(default [1 0.16 0]) 018 % 019 % This is a transformation model where the transformation is chosen to 020 % be a monotonic exponential function, calibrated such that the first 021 % 3 moments of the transformed model G(y)=g^-1(y) match the moments of 022 % the true process. However, the skewness is limited by ABS(SKEW)<2.82. 023 % Information about the moments of the process can be 024 % obtained by site specific data, laboratory measurements or by resort 025 % to theoretical models (see spec2skew). According to Ochi it is 026 % appropriate for a process with very strong non-linear characteristics. 027 % 028 % g(x) = ((1-exp(-gamma*(x-mean)/sigma))/gamma-mean2)/sigma2 029 % where 030 % gamma = 1.28*a for x>=mean 031 % 3*a otherwise 032 % mean, 033 % sigma = standard deviation and mean, respectively, of the process. 034 % mean2, 035 % sigma2 = normalizing parameters in the transformed world, i.e., to 036 % make the gaussian process in the transformed world is 037 % N(0,1). 038 % 039 % The unknown parameters a, mean2 and sigma2 are found by solving the 040 % following non-linear equations: 041 % 042 % a*(sigma2^2+mean2^2)+mean2 = 0 043 % sigma2^2-2*a^2*sigma2^4 = 1 044 % 2*a*sigma2^4*(3-8*a^2*sigma2^2) = skew 045 % 046 % NOTE: - by specifying NaN's in the data vector default values will be used. 047 % - if length(data) is shorter than the parameters needed then the 048 % default values are used for the parameters not specified. 049 % - The gaussian process in the transformed world is N(0,1) 050 % - g does not have continous derivatives of 2'nd order or higher. 051 % 052 % Example: Simulate a Transformed Gaussian process: 053 % Hm0=7;Tp=11; 054 % S = jonswap([],[Hm0 Tp]); [sk ku ma]=spec2skew(S); 055 % g = ochitr([],[Hm0/4,sk,ma]); g2=[g(:,1), g(:,2)*Hm0/4]; 056 % ys = spec2sdat(S,15000); % Simulated in the Gaussian world 057 % xs = gaus2dat(ys,g2); % Transformed to the real world 058 % 059 % See also spec2skew, hermitetr, lc2tr, dat2tr 060 061 062 % References: 063 % Ochi, M.K. and Ahn, K. (1994) 064 % 'Non-Gaussian probability distribution of coastal waves.' 065 % In Proc. 24th Conf. Coastal Engng, Vol. 1, pp 482-496 066 % 067 % Michel K. Ochi (1998), 068 % "OCEAN WAVES, The stochastic approach", 069 % OCEAN TECHNOLOGY series 6, Cambridge, pp 255-275. 070 071 072 % tested on matlab 5.1 073 % History: 074 % revised pab 02.01.2001 075 % - default x is now levels([-5 5 513])*sa+ma -> better to have the discretization 076 % represented with exact numbers, especially when calculating 077 % derivatives of the transformation numerically. 078 % - added fmins, fzero 079 % - moved some code into ochifun. 080 % revised pab 24.05.2000 081 % - removed inline object with string object 082 % by pab 03.03.2000 083 084 data2=[1 0.16 0]; % default values 085 if nargin>=2 & any(~isnan(data)) 086 ind=find(~isnan(data(1:min(length(data),3)))); 087 data2(ind)=data(ind); 088 end 089 sig1=data2(1); skew=data2(2); ma=data2(3); 090 091 if abs(skew)>2.82842712474619, 092 error('Skewness must be less than 2.82842') 093 end 094 095 if nargin<1|isempty(y), y=linspace(-5*sig1+ma,5*sig1+ma,513); end 096 097 %disp('1') 098 phat = ochitrfit(sig1,skew,ma); 099 %disp('2') 100 101 102 %disp('3') 103 if nargout>1, 104 [g,t0] = ochifun(y,phat); 105 else 106 g = ochifun(y,phat); 107 end 108 return 109 110 111 112 113 function phat=ochitrfit(sig1,skew,ma) 114 115 if skew==0 116 phat = [0, 0, sig1, ma,1 0]; 117 return 118 end 119 % Solve the equations to obtain the gamma parameters: 120 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 121 % a*(sig2^2+ma2^2)+ma2 = 0 122 % sig2^2-2*a^2*sig2^4 = E(y^2) % =1 123 % 2*a*sig2^4*(3-8*a^2*sig2^2) = E(y^3) % = skew 124 125 % Let x = [a sig2^2 ] 126 % Set up the 2D non-linear equations for a and sig2^2: 127 eqstr='[x(2)-2.*x(1).^2.*x(2).^2-P1, 2.*x(1).*x(2).^2.*(3-8.*x(1).^2.*x(2))-P2 ]'; 128 % Or solve the following 1D non-linear equation for sig2^2: 129 eqstr2 = '-sqrt(abs(x-P1)*2).*(3.*x-4*abs(x-P1))+abs(P2)'; 130 %g3 = inline(eqstr2,2); 131 132 g2 = eqstr2; 133 g1 = eqstr; 134 135 v = version; ix = find(v=='.'); 136 vnr= str2num(v(1:ix(min(2,length(ix)))-1)); 137 138 % start value for: a sig2^2 139 X0=[0.01 sig1^2]; 140 Xi =[1 2]; % Start interval where sig2^2 is located. 141 if vnr>5.2 142 % sol = fsolve(g1,X0,[],1,skew);%sig1.^2,skew*sig1^3); 143 opt = [];%optimset('disp','iter'); 144 if 1, 145 sig22 = fzero(g2,Xi,opt,1,skew); % smallest solution for sig22 146 a = sign(skew)*sqrt(abs(sig22-1)/2/sig22^2); 147 sol = [a sig22]; 148 else 149 % find the solution by least squares 150 g1 = [ 'sum(' g1 '.^2)']; 151 sol = fminsearch(g1,X0,opt,1,skew);%sig1.^2,skew*sig1^3); 152 end 153 else 154 % sol = fsolve(g1,X0,[],[],1,skew);%sig1.^2,skew*sig1^3); 155 trace1=[]; 156 if 1, 157 sig22 = fzero(g2,[1 2],[],trace1,1,skew); 158 a = sign(skew)*sqrt(abs(sig22-1)/2/sig22^2); 159 sol = [a sig22]; 160 else 161 % find the solution by least squares 162 g1 = [ 'sum(' g1 '.^2)']; 163 %sol = fmins(g1,X0,[],1,sig1.^2,skew*sig1^3); 164 sol = fmins(g1,X0,[],trace1,1.^2,skew*1^3); 165 end 166 end 167 168 a = sol(1); 169 sig22 = sol(2); 170 sig2 = sqrt(sig22); 171 172 % Solve the following 2nd order equation to obtain ma2 173 % a*(sig2^2+ma2^2)+ma2 = 0 174 my2 = (-1-sqrt(1-4*a^2*sig22))./a; % Largest mean 175 ma2 = a*sig22/my2 ; % choose the smallest mean 176 gam_ab = [1.28 3]*a; % this is valid for processes with very strong 177 % nonlinear characteristics 178 phat = [gam_ab sig1 ma sig2 ma2]; 179 180 181 return 182 183 184 185 186 187 188
Comments or corrections to the WAFO group