OCHIFUN Calculates the 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] = ochifun(y,data); 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 of the data from the Gaussian model. x = a row vector with x-values. (default linspace(-5*sigma,5*sigma,513)+mean) data = [gamma_a gamma_b sigma mean sigma2 mean2], transformation parameters, standard deviation, and mean, respectively. (default [0.1 0.15 1 0 1 0]) Mean2 and sigma2 are normalizing parameters in the transformed world, i.e., to make the gaussian process in the transformed world is N(0,1). This is a transformation model where the transformation is chosen to be a monotonic exponential function: g(x) = ((1-exp(-gamma*(x-mean)/sigma))/gamma-mean2)/sigma2 where gamma = gamma_a for x>=mean gamma_b otherwise According to Ochi it is appropriate for a process with very strong non-linear characteristics. NOTE: - g does not have continous derivatives of 2'nd order or higher unless gamma_a==gamma_b. See also ochitr, hermitetr, lc2tr, dat2tr
Calculates transformation, g, proposed by Ochi et al. |
001 function [g,t0]=ochifun(y,data) 002 %OCHIFUN Calculates the transformation g proposed by Ochi et al. 003 % Assumption: a Gaussian process, Y, is related to the 004 % non-Gaussian process, X, by Y = g(X). 005 % 006 % CALL: [g, test] = ochifun(y,data); 007 % 008 % g = [x g(x)], a two column matrix with the transformation g(x). 009 % test = int (g(x)-x)^2 dx where int. limits is given by X. This 010 % is a measure of departure of the data from the Gaussian model. 011 % x = a row vector with x-values. 012 % (default linspace(-5*sigma,5*sigma,513)+mean) 013 % data = [gamma_a gamma_b sigma mean sigma2 mean2], 014 % transformation parameters, standard deviation, and mean, 015 % respectively. (default [0.1 0.15 1 0 1 0]) 016 % Mean2 and sigma2 are normalizing parameters in the 017 % transformed world, i.e., to make the gaussian process in 018 % the transformed world is N(0,1). 019 % 020 % This is a transformation model where the transformation is chosen to 021 % be a monotonic exponential function: 022 % 023 % g(x) = ((1-exp(-gamma*(x-mean)/sigma))/gamma-mean2)/sigma2 024 % where 025 % gamma = gamma_a for x>=mean 026 % gamma_b otherwise 027 % 028 % According to Ochi it is appropriate for a process with very strong 029 % non-linear characteristics. 030 % 031 % NOTE: - g does not have continous derivatives of 2'nd order or higher 032 % unless gamma_a==gamma_b. 033 % 034 % See also ochitr, hermitetr, lc2tr, dat2tr 035 036 037 % References: 038 % Ochi, M.K. and Ahn, K. (1994) 039 % 'Non-Gaussian probability distribution of coastal waves.' 040 % In Proc. 24th Conf. Coastal Engng, Vol. 1, pp 482-496 041 % 042 % Michel K. Ochi (1998), 043 % "OCEAN WAVES, The stochastic approach", 044 % OCEAN TECHNOLOGY series 6, Cambridge, pp 255-275. 045 046 047 % tested on matlab 5.1 048 % History: 049 % revised pab 02.01.2001 050 % - changed name to ochifun 051 % - fixed a bug: gamma =0 is now handled correctly 052 % - default x is now levels([-5 5 513])*sa+ma -> better to have the 053 % discretization 054 % represented with exact numbers, especially when calculating 055 % derivatives of the transformation numerically. 056 % revised pab 21.02.2000 057 % - added ma ma2 058 % - changed name to ochitr 059 % - added references 060 % - changed default value for y 061 % - fixed a normalization bug 062 063 %1./data2(1:2) 064 065 data2 = [0.1 0.15 1 0 1 0]; 066 if nargin>=2 & any(~isnan(data)) 067 ind=find(~isnan(data(1:min(length(data),6)))); 068 data2(ind)=data(ind); 069 end 070 ga =data2(1); gb = data2(2); 071 sigma = data2(3); ma = data2(4); 072 sigma2 = data2(5); ma2 = data2(6); 073 if nargin<1|isempty(y); 074 y = linspace(-5*sigma+ma,5*sigma+ma,513)'; 075 else 076 y=y(:); 077 end 078 g=zeros(length(y),2); 079 g(:,1)=y; 080 igp=(y>=ma); 081 igm=find(~igp); 082 083 yn = (y-ma)/sigma; 084 if ga==0, 085 g(igp,2)=yn(igp); 086 else 087 g(igp,2)=(1-exp(-ga*yn(igp)))/ga; 088 end 089 if gb==0, 090 g(igm,2)=yn(igm); 091 else 092 g(igm,2)=(1-exp(-gb*yn(igm)))/gb; 093 end 094 095 g(:,2)=(g(:,2)-ma2)/sigma2; 096 097 if nargout>1 098 t0 = trapz(yn,(yn-g(:,2)).^2); 099 end 100
Comments or corrections to the WAFO group