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
Simulates a Gaussian process and its derivative | |
Estimate transformation, g, from data. | |
Computes covariance function and its derivatives | |
Create or alter TRANSFORM OPTIONS structure. | |
Display message and abort function. | |
Hold current graph. | |
Convert number to string. (Fast version) | |
Linear plot. | |
X-axis label. | |
Y-axis label. |
% CHAPTER2 Modelling random loads and stochastic waves | |
Variability of simulated e(g(u)-u) (circles) |
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
Comments or corrections to the WAFO group