Home > wafo > trgauss > mctrtest.m

mctrtest

PURPOSE ^

Test if a stochastic process is Gaussian.

SYNOPSIS ^

test2 = mctrtest(S,Np,test0,def,opt)

DESCRIPTION ^

 MCTRTEST Test if a stochastic process is Gaussian.
 
  CALL:  test1 = mctrtest(S,[Np,Ns],test0,def,options);
 
  test1,
     test0 = simulated and observed value of e(g)=int (g(u)-u)^2 du,
             respectively, where int limits is given by OPTIONS.PARAM. 
          
       S   = spectral density structure 
 
        Np = # of points simulated
        Ns = # of independent simulations (default  100)
 
    def    = 'nonlinear' : transform based on smoothed crossing intensity (default)
             'mnonlinear': transform based on smoothed marginal distribution
   options = options structure defining how the estimation of the
             transformation is done. (default troptset('dat2tr'))
 
  MCTRTEST simulates  e(g(u)-u) = int (g(u)-u)^2 du  for Gaussian processes 
  given the spectral density, S. The result is plotted if test0 is given.
  This is useful for testing if the process X(t) is Gaussian.
  If 95% of TEST1 is less than TEST0 then X(t) is not Gaussian at a 5% level.
 
  Example:
  Hm0 = 7;
  S0 = jonswap([],Hm0); g=ochitr([],[Hm0/4]); S=S0;
  S.tr=g;S.tr(:,2)=g(:,2)*Hm0/4;
  xs = spec2sdat(S,2^13);
  [g0 t0] = dat2tr(xs);
  t1 = mctrtest(S0,[2^13 50],t0); 
 
  See also  cov2sdat, dat2tr, troptset

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function test2 = mctrtest(S,Np,test0,def,opt)
002 %MCTRTEST Test if a stochastic process is Gaussian.
003 %
004 % CALL:  test1 = mctrtest(S,[Np,Ns],test0,def,options);
005 %
006 % test1,
007 %    test0 = simulated and observed value of e(g)=int (g(u)-u)^2 du,
008 %            respectively, where int limits is given by OPTIONS.PARAM. 
009 %         
010 %      S   = spectral density structure 
011 %
012 %       Np = # of points simulated
013 %       Ns = # of independent simulations (default  100)
014 %
015 %   def    = 'nonlinear' : transform based on smoothed crossing intensity (default)
016 %            'mnonlinear': transform based on smoothed marginal distribution
017 %  options = options structure defining how the estimation of the
018 %            transformation is done. (default troptset('dat2tr'))
019 %
020 % MCTRTEST simulates  e(g(u)-u) = int (g(u)-u)^2 du  for Gaussian processes 
021 % given the spectral density, S. The result is plotted if test0 is given.
022 % This is useful for testing if the process X(t) is Gaussian.
023 % If 95% of TEST1 is less than TEST0 then X(t) is not Gaussian at a 5% level.
024 %
025 % Example:
026 % Hm0 = 7;
027 % S0 = jonswap([],Hm0); g=ochitr([],[Hm0/4]); S=S0;
028 % S.tr=g;S.tr(:,2)=g(:,2)*Hm0/4;
029 % xs = spec2sdat(S,2^13);
030 % [g0 t0] = dat2tr(xs);
031 % t1 = mctrtest(S0,[2^13 50],t0); 
032 %
033 % See also  cov2sdat, dat2tr, troptset
034 
035 %Tested on: Matlab 5.3, 5.2, 5.1
036 % History:
037 % revised pab jan2005
038 % changed order of input so that chapter1.m works
039 % revised pab Feb2004
040 % -changed h1line  
041 % revised pab 29.12.2000
042 % - added options and def to input due to changed calling syntax for dat2tr.
043 % - updated the help header
044 % revised by pab 12.11.99
045 %    fixed a bug, string input to dat2tr 'nonlinear' 
046 % revised by pab 12.10.99
047 %    updated help header 
048 % revised by pab 11.08.99
049 % changed name from mctest to mctrtest
050 % by pab 11.11.98
051 
052 maxsize = 200000; % must divide the computations due to limited memory
053 
054 
055 error(nargchk(2,5,nargin))
056 if nargin<4|isempty(def),
057   def = 'nonlinear';
058 end
059 if nargin<5|isempty(opt),
060   opt = troptset('dat2tr');
061 end
062 opt = troptset(opt,'multip',1);
063 
064 if nargin<3|isempty(test0)
065   plotflag=0;
066 else 
067   plotflag=1;
068 end
069 
070 Ns=100;
071 nnp=length(Np);
072 if nnp>=2 
073   Ns=Np(2);
074 end
075 Np=Np(1);
076 if Ns>50
077   disp('  ... be patient this may take a while')
078 end
079 
080 test1 = [];
081 rep = floor(Np*Ns/maxsize)+1;
082 
083 Nstep = floor(Ns/rep);
084 
085 R = spec2cov(S);
086 
087 for ix=1:rep,
088   xs = cov2sdat(R,[Np Nstep]);
089   [g, tmp] = dat2tr(xs,def,opt);
090   test1 = [test1; tmp(:)];
091   disp(['finished ' num2str(ix) ' of ' num2str(rep)] )
092 end
093 if rep>1,
094   xs = cov2sdat(R,[Np rem(Ns,rep)]);
095   [g tmp] = dat2tr(xs,def,opt);
096   test1 = [test1; tmp(:)];
097 end
098 
099 if (nargout>0 | plotflag==0),
100   test2=test1;
101 end
102 
103 
104 if plotflag 
105   plot(test1,'o'),hold on
106   if 1 
107     plot([1 Ns],test0*[1 1],'--'),
108   end
109   hold off
110   ylabel('e(g(u)-u)')
111   xlabel('Simulation number')
112 end
113

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