Home > wafo > trgauss > ochitr.m

ochitr

PURPOSE ^

Calculates transformation, g, proposed by Ochi et al.

SYNOPSIS ^

[g,t0]=ochitr(y,data)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003