Home > wafo > trgauss > cdf2tr.m

cdf2tr

PURPOSE ^

Estimate transformation, g, from observed CDF.

SYNOPSIS ^

[g, test, g2] = cdf2tr(Fx1,ma ,sa,varargin)

DESCRIPTION ^

 CDF2TR Estimate transformation, g, from observed CDF. 
  
         Assumption: a Gaussian process, Y, is related to the 
                     non-Gaussian process, X, by Y = g(X).  
   
   CALL [g,test,g2] = cdf2tr(F,ma,sa,options); 
  
      g,g2  = smoothed and empirical estimate of the transformation  g.      
      test  = test observator int (g(u)-u)^2 du  where int limits is 
              given by OPTIONS.PARAM. This is a measure of departure of the  
              data from the Gaussian model. 
      F     = empirical CDF of X(t), a 2 column matrix. 
      ma,sa = mean and standard deviation of the process X(t). 
    options = options structure defining how the smoothing is done. 
              (See troptset for default values) 
  
     The empirical CDF is usually very irregular. 
   More than one local maximum of the empirical CDF 
   may cause poor fit of the transformation. In such case one 
   should use a smaller value of GSM or set a larger variance for GVAR.  
     If X(t) is likely to cross levels higher than 5 standard deviations   
   then the vector param has to be modified.  For example if X(t) is  
   unlikely to cross a level of 7 standard deviations one can use  
   param = [-7 7 513]. 
  
  Example 
  Hm0 = 7; 
  S = jonswap([],Hm0); g=ochitr([],[Hm0/4]);  
  S.tr = g; S.tr(:,2)=g(:,2)*Hm0/4; 
  xs = spec2sdat(S,2^13); 
  Fx = empdistr(xs(:,2)); 
  g0 = cdf2tr(Fx,0,Hm0/4,troptset('plot',1));   % Plot final estimate 
  g1 = lc2tr(Fx,0,Hm0/4,troptset('gvar', .5 )); % More weight on all points 
  g2 = lc2tr(Fx,0,Hm0/4,'gvar', [3.5 .5 3.5]);  % Less weight on the ends 
  hold on, trplot(g1,g)                                   % Check the fit 
  trplot(g2) 
  
  See also  troptset, lc2tr

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [g, test, g2] = cdf2tr(Fx1,ma ,sa,varargin) 
002 %CDF2TR Estimate transformation, g, from observed CDF. 
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,g2] = cdf2tr(F,ma,sa,options); 
008 % 
009 %     g,g2  = smoothed and empirical estimate of the transformation  g.      
010 %     test  = test observator int (g(u)-u)^2 du  where int limits is 
011 %             given by OPTIONS.PARAM. This is a measure of departure of the  
012 %             data from the Gaussian model. 
013 %     F     = empirical CDF of X(t), a 2 column matrix. 
014 %     ma,sa = mean and standard deviation of the process X(t). 
015 %   options = options structure defining how the smoothing is done. 
016 %             (See troptset for default values) 
017 % 
018 %    The empirical CDF is usually very irregular. 
019 %  More than one local maximum of the empirical CDF 
020 %  may cause poor fit of the transformation. In such case one 
021 %  should use a smaller value of GSM or set a larger variance for GVAR.  
022 %    If X(t) is likely to cross levels higher than 5 standard deviations   
023 %  then the vector param has to be modified.  For example if X(t) is  
024 %  unlikely to cross a level of 7 standard deviations one can use  
025 %  param = [-7 7 513]. 
026 % 
027 % Example 
028 % Hm0 = 7; 
029 % S = jonswap([],Hm0); g=ochitr([],[Hm0/4]);  
030 % S.tr = g; S.tr(:,2)=g(:,2)*Hm0/4; 
031 % xs = spec2sdat(S,2^13); 
032 % Fx = empdistr(xs(:,2)); 
033 % g0 = cdf2tr(Fx,0,Hm0/4,troptset('plot',1));   % Plot final estimate 
034 % g1 = lc2tr(Fx,0,Hm0/4,troptset('gvar', .5 )); % More weight on all points 
035 % g2 = lc2tr(Fx,0,Hm0/4,'gvar', [3.5 .5 3.5]);  % Less weight on the ends 
036 % hold on, trplot(g1,g)                                   % Check the fit 
037 % trplot(g2) 
038 % 
039 % See also  troptset, lc2tr 
040  
041 % History 
042 % revised Feb2004   
043 % Revised pab Dec2003 
044 %  fixed a bug F -> Fx1 
045 % by pab 29.12.2000 
046 % - default param is now [-5 5 513] -> better to have the discretization 
047 %  represented with exact numbers, especially when calculating 
048 %  derivatives of the transformation numerically. 
049  
050 opt = troptset('chkder','on','plotflag','off','gsm',.05,.... 
051     'param',[-5 5 513],'delay',2,'linextrap','on','ntr',1000,'ne',7,'gvar',1); 
052 % If just 'defaults' passed in, return the default options in g 
053 if nargin==1 & nargout <= 1 & isequal(Fx1,'defaults') 
054   g = opt;  
055   return 
056 end 
057 error(nargchk(3,inf,nargin)) 
058 if nargin>=4,  opt=troptset(opt,varargin{:}); end 
059 switch opt.chkder; 
060   case 'off', chkder = 0; 
061   case 'on',  chkder = 1; 
062   otherwise,  chkder = opt.chkder; 
063 end 
064 switch opt.linextrap; 
065   case 'off', def = 0; 
066   case 'on',  def = 1; 
067   otherwise,  def = opt.linextrap; 
068 end 
069  
070 switch opt.plotflag 
071   case {'none','off'},   plotflag = 0; 
072   case 'final', plotflag = 1; 
073   case 'iter',  plotflag = 2; 
074   otherwise,    plotflag = opt.plotflag; 
075 end 
076  
077  
078 Ne = opt.ne; 
079 if length(Fx1)>opt.ntr & opt.ntr>0 
080   x0 = linspace(Fx1(1+Ne,1),Fx1(end-Ne,1),opt.ntr)'; 
081   Fx = [ x0,interp1q(Fx1(:,1),Fx1(:,2),x0)]; 
082   Ne=0; 
083 else 
084   Fx = Fx1; 
085 end 
086 uu = levels(opt.param)'; 
087 g = [sa*uu+ma zeros(opt.param(3),1)]; 
088 ncr = length(Fx); 
089  
090  
091 ng = length(opt.gvar); 
092 if ng==1 
093   gvar = opt.gvar(ones(ncr,1)); 
094 else 
095   gvar = interp1(linspace(0,1,ng)',opt.gvar(:),linspace(0,1,ncr)','*linear');   
096 end 
097   
098  
099 ind = find(diff(Fx(:,1))>0); % remove equal points 
100 ind1 = ind(Ne+1:end-Ne);   
101 tmp = wnorminv(Fx(ind,2)); 
102  
103 g(:,2) = smooth(Fx(ind1,1),tmp(Ne+1:end-Ne),opt.gsm,g(:,1),def,gvar); 
104  
105 if chkder~=0 
106   for ix = 1:5 
107     dy = diff(g(:,2)); 
108     if any(dy<=0) 
109       warning('The empirical distribution is not sufficiently smoothed.') 
110       disp('        The estimated transfer function, g, is not ') 
111       disp('        a strictly increasing function.') 
112       dy(dy>0)=eps; 
113       gvar = -([dy;0]+[0;dy])/2+eps; 
114       g(:,2) = smooth(g(:,1),g(:,2),1,g(:,1),def,ix*gvar); 
115     else  
116       break 
117     end 
118   end 
119 end 
120 if nargout>1 
121   test = sqrt(trapz(uu,(uu-g(:,2)).^2));   
122 end 
123 if nargout>2 
124   g2 = [Fx(ind,1) tmp]; 
125 end 
126 if plotflag>0 
127   trplot(g,g2,ma,sa) 
128   if plotflag>1 
129     pause(opt.delay) 
130   end 
131 end 
132  
133  
134  
135

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