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)

## CROSS-REFERENCE INFORMATION

This function calls:
 levels Calculates discrete levels given the parameter matrix. smooth Calculates a smoothing spline. troptset Create or alter TRANSFORM OPTIONS structure. trplot Plots transformation, g, eg. estimated with dat2tr. wnorminv Inverse of the Normal distribution function diff Difference and approximate derivative. error Display message and abort function. interp1 1-D interpolation (table lookup) interp1q Quick 1-D linear interpolation. isequal True if arrays are numerically equal. linspace Linearly spaced vector. pause Wait for user response. trapz Trapezoidal numerical integration. warning Display warning message; disable or enable warning messages.
This function is called by:
 dat2tr Estimate transformation, g, from data.

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