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);

## CROSS-REFERENCE INFORMATION

This function calls:
 cov2sdat Simulates a Gaussian process and its derivative dat2tr Estimate transformation, g, from data. spec2cov Computes covariance function and its derivatives troptset Create or alter TRANSFORM OPTIONS structure. error Display message and abort function. hold Hold current graph. num2str Convert number to string. (Fast version) plot Linear plot. xlabel X-axis label. ylabel Y-axis label.
This function is called by:
 Chapter2 % CHAPTER2 Modelling random loads and stochastic waves recfig6 Variability of simulated e(g(u)-u) (circles)

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