Home > wafo > trgauss > hermitetr.m

hermitetr

PURPOSE ^

Calculate transformation, g, proposed by Winterstein

SYNOPSIS ^

[g ,t0]=hermitetr(x,data,def)

DESCRIPTION ^

 HERMITETR Calculate transformation, g, proposed by Winterstein 
  
            Assumption: a Gaussian process, Y, is related to the 
                      non-Gaussian process, X, by Y = g(X).  
  
   CALL:  [g,test] = hermitetr(x,data,def); 
          [g,test] = hermitetr(x,S,def); 
  
         g    = [x g(x)] a two column matrix with the transformation g(x). 
         test = int (g(x)-x)^2 dx  where int. limits is given by X. This 
                is a measure of departure from the Gaussian model. 
         x    = a row vector with x-values.  
                (default linspace(-5*sigma,5*sigma,501)+mean) 
         data = [sigma skew kurt mean] is the  standard deviation, 
                skewness, kurtosis and mean of the process, 
                respectively. skew=kurt-3=0 for a Gaussian process. 
                This is fairly accurate if kurt>=0 and  
                0<=skew^2 <= 8*kurt/9  (default  [1 0.16 3.04 0]) 
         S    = spectral density struct from which  
                [sigma skew kurt mean] is calculated using spec2skew. 
         def  = 1  Winterstein et. al. (1994) parametrization (default) 
                2  Winterstein (1988) parametrization 
  
  HERMITETR is a hermite transformation model where the transformation is 
  chosen to be monotonic cubic polynomial, calibrated such that the first  
  4 moments of the transformed model G(y)=g^-1(y) match the moments of 
  the true process. Information about the moments of the process can be 
  obtained by site specific data, laboratory measurements or by resort to 
  theoretical models (see spec2skew). 
  
  If kurt<3 (hardening model) 
     g(x) =  xn - c3(xn^2-1) - c4*(xn^3-3*xn)  
  where  
    xn = (x-mean)/sigma 
    c3 = skew/6 
    c4 = (kurt-3)/24 
  
  If kurt>=3 (softening model) 
     G(y) = mean + K*sigma*[ y + c3(y^2-1) + c4*(y^3-3*y) ] 
   where 
     y  = g(x) = G^-1(x) 
     K  = 1/sqrt(1+2*c3^2+6*c4^2) 
  If def = 2 : 
     c3 = skew/(6*(1+6*c4)) 
     c4 = [sqrt(1+1.5*(kurt-3))-1]/18  
  If def = 1 : 
     c3  = skew/6*(1-0.015*abs(skew)+0.3*skew^2)/(1+0.2*(kurt-3)) 
     c4  = 0.1*((1+1.25*(kurt-3))^(1/3)-1)*c41 
     c41 = (1-1.43*skew^2/(kurt-3))^(1-0.1*(kurt)^0.8) 
  
   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) 
  
  Example: Simulate a Transformed Gaussian process: 
   Hm0=7;Tp=11; 
   S = jonswap([],[Hm0 Tp]); g=hermitetr*Hm0/4;  
   ys = spec2sdat(S,15000);   % Simulated in the Gaussian world 
   xs = gaus2dat(ys,g);      % Transformed to the real world 
  
  See also  spec2skew, ochitr, lc2tr, dat2tr

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [g ,t0]=hermitetr(x,data,def) 
002 %HERMITETR Calculate transformation, g, proposed by Winterstein 
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] = hermitetr(x,data,def); 
008 %         [g,test] = hermitetr(x,S,def); 
009 % 
010 %        g    = [x g(x)] a two column matrix with the transformation g(x). 
011 %        test = int (g(x)-x)^2 dx  where int. limits is given by X. This 
012 %               is a measure of departure from the Gaussian model. 
013 %        x    = a row vector with x-values.  
014 %               (default linspace(-5*sigma,5*sigma,501)+mean) 
015 %        data = [sigma skew kurt mean] is the  standard deviation, 
016 %               skewness, kurtosis and mean of the process, 
017 %               respectively. skew=kurt-3=0 for a Gaussian process. 
018 %               This is fairly accurate if kurt>=0 and  
019 %               0<=skew^2 <= 8*kurt/9  (default  [1 0.16 3.04 0]) 
020 %        S    = spectral density struct from which  
021 %               [sigma skew kurt mean] is calculated using spec2skew. 
022 %        def  = 1  Winterstein et. al. (1994) parametrization (default) 
023 %               2  Winterstein (1988) parametrization 
024 % 
025 % HERMITETR is a hermite transformation model where the transformation is 
026 % chosen to be monotonic cubic polynomial, calibrated such that the first  
027 % 4 moments of the transformed model G(y)=g^-1(y) match the moments of 
028 % the true process. Information about the moments of the process can be 
029 % obtained by site specific data, laboratory measurements or by resort to 
030 % theoretical models (see spec2skew). 
031 % 
032 % If kurt<3 (hardening model) 
033 %    g(x) =  xn - c3(xn^2-1) - c4*(xn^3-3*xn)  
034 % where  
035 %   xn = (x-mean)/sigma 
036 %   c3 = skew/6 
037 %   c4 = (kurt-3)/24 
038 % 
039 % If kurt>=3 (softening model) 
040 %    G(y) = mean + K*sigma*[ y + c3(y^2-1) + c4*(y^3-3*y) ] 
041 %  where 
042 %    y  = g(x) = G^-1(x) 
043 %    K  = 1/sqrt(1+2*c3^2+6*c4^2) 
044 % If def = 2 : 
045 %    c3 = skew/(6*(1+6*c4)) 
046 %    c4 = [sqrt(1+1.5*(kurt-3))-1]/18  
047 % If def = 1 : 
048 %    c3  = skew/6*(1-0.015*abs(skew)+0.3*skew^2)/(1+0.2*(kurt-3)) 
049 %    c4  = 0.1*((1+1.25*(kurt-3))^(1/3)-1)*c41 
050 %    c41 = (1-1.43*skew^2/(kurt-3))^(1-0.1*(kurt)^0.8) 
051 % 
052 %  NOTE: - by specifying NaN's in the data vector default values will be used. 
053 %        - if length(data) is shorter than the parameters needed then the 
054 %          default values are used for the parameters not specified.  
055 %        - The gaussian process in the transformed world is N(0,1) 
056 % 
057 % Example: Simulate a Transformed Gaussian process: 
058 %  Hm0=7;Tp=11; 
059 %  S = jonswap([],[Hm0 Tp]); g=hermitetr*Hm0/4;  
060 %  ys = spec2sdat(S,15000);   % Simulated in the Gaussian world 
061 %  xs = gaus2dat(ys,g);      % Transformed to the real world 
062 % 
063 % See also  spec2skew, ochitr, lc2tr, dat2tr 
064  
065 % References: 
066 % Winterstein, S.R. (1988) 
067 % 'Nonlinear vibration models for extremes and fatigue.' 
068 % J. Engng. Mech., ASCE, Vol 114, No 10, pp 1772-1790 
069 % 
070 % Winterstein, S.R, Ude, T.C. and Kleiven, G. (1994) 
071 % "Springing and slow drift responses: 
072 %  predicted extremes and fatigue vs. simulation" 
073 % In Proc. 7th International behaviour of Offshore structures, (BOSS) 
074 % Vol. 3, pp.1-15 
075  
076  
077  
078 % tested on matlab 5.2 
079 % History: 
080 % revised pab dec 2003   
081 % revised pab 02.04.2001 
082 % -Added Spectrum as a possible input 
083 %revised pab 02.01.2000 
084 % - added t0, check on that g is monotone 
085 % - moved some code to hermitefun 
086 % by pab 21.02.2000 
087  
088  
089 data2=[1 0.16 3.04 0]; % default values 
090 if nargin>=2 & ~isempty(data) 
091   switch class(data) 
092     case 'double', 
093      if  any(~isnan(data)) 
094        ind=find(~isnan(data(1:min(length(data),4)))); 
095        data2(ind)=data(ind); 
096      end 
097    case 'struct', % data is a spctral density struct 
098     [skew, kurt, ma, sigma] = spec2skew(data); 
099     data2 =  [sigma skew, kurt, ma]; 
100    otherwise 
101     warning('Wrong input data') 
102     end 
103 end 
104 sigma=data2(1); skew=data2(2); kurt=data2(3);  ma=data2(4); 
105 if nargin<1|isempty(x),   x   = linspace(-5*sigma,5*sigma,501)+ma; end 
106 if nargin<3|isempty(def), def = 1;end 
107  
108 %skew,ga2 
109 ga2 = kurt-3; 
110 if ga2<0 
111   c4 = ga2/24; 
112   c3 = skew/6; 
113 else 
114   switch def 
115     case 2, % Winterstein 1988 parametrization 
116       if skew^2>8*(ga2+3)/9, 
117     disp('warning: kurtosis too low compared to the skewness') 
118       end 
119       c4 = (sqrt(1+1.5*ga2)-1)/18;  
120       c3 = skew/(6*(1+6*c4)); 
121     otherwise,   % Winterstein et. al. 1994 parametrization intended to 
122       % apply for the range:  0 <= ga2 < 12 and 0<= skew^2 < 2*ga2/3 
123       if skew^2>2*(ga2)/3, 
124     disp('Warning: kurtosis too low compared to the skewness') 
125       end 
126       if (ga2 < 0)| (12 < ga2) 
127     disp('Warning: kurtosis must be between 0 and 12') 
128       end 
129       c3 = skew/6*(1-0.015*abs(skew)+0.3*skew^2)/(1+0.2*ga2); 
130       if ga2==0, 
131     c4=0; 
132       else 
133     c41= (1-1.43*skew.^2./ga2).^(1-0.1*(ga2+3).^0.8); 
134     c4 = 0.1*((1+1.25*ga2)^(1/3)-1)*c41; 
135       end 
136   end 
137 end 
138 if ~isreal(c4)|~isreal(c3) 
139  error('Unable to calculate the polynomial') 
140 end 
141  
142 if nargout>1 
143   [g, t0]=hermitefun([c3 c4],x(:),ma,sigma,ga2); 
144 else 
145   g = hermitefun([c3 c4],x(:),ma,sigma,ga2); 
146 end 
147 return 
148  
149  
150  
151  
152  
153

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