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

## CROSS-REFERENCE INFORMATION

This function calls:
 ochifun Calculates the transformation g proposed by Ochi et al. error Display message and abort function. fmins fminsearch Multidimensional unconstrained nonlinear minimization (Nelder-Mead). fzero Scalar nonlinear zero finding. linspace Linearly spaced vector. str2num Convert string matrix to numeric array. version MATLAB version number.
This function is called by:
 dat2tr Estimate transformation, g, from data.

## 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 %
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