line source
+ − % SFtest perfomes a Spectral F-Test on PSDs.
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ − %
+ − % DESCRIPTION: SFtest performes a Spectral F-Test on two input PSD objects.
+ − % The null hypothesis H0 (the two PSDs belong to the same
+ − % statistical distribution) is rejected at the confidence
+ − % level for the alternative hypotheis H1 (the two PSDs belong
+ − % to different statistical distributions) if the test
+ − % statistic falls in the critical region.
+ − % SFtest uses utils.math.Ftest which does the test at each
+ − % frequency bin.
+ − %
+ − % VERSION: $Id: SFtest.m,v 1.3 2011/03/07 17:38:05 congedo Exp $
+ − %
+ − % HISTORY: 18-02-2011 G. Congedo
+ − %
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ −
+ − function test = SFtest(X,Y,alpha,showPlots)
+ −
+ − % Assume a two-tailed test
+ − twoTailed = 1;
+ −
+ − % Extract the degree of freedom
+ − dofX = X.getdof.y;
+ − dofY = Y.getdof.y;
+ −
+ − % Ratio of the two spectra: this test statistic is F-distributed
+ − F = X/Y;
+ − F.setYunits('');
+ − F.setName('F statistic');
+ −
+ − n = numel(F.y);
+ −
+ − % Interquartile range
+ − Fy = sort(F.y);
+ − m = median(Fy);
+ − q1 = median(Fy(Fy<=m));
+ − q3 = median(Fy(Fy>=m));
+ − iqr = q3 - q1;
+ −
+ − % Freedman–Diaconis rule
+ − binSz = 2*iqr*n^(-1/3);
+ − nBin = round((max(Fy)-min(Fy))/binSz);
+ −
+ − % Sample PDF
+ − samplePDF = hist(F,plist('N',nBin,'norm',1));
+ − samplePDF.setName('sample PDF');
+ −
+ − % Sample PDF moments
+ − % sampleK = utils.math.Kurt(samplePDF.x);
+ − % sampleS = utils.math.Skew(samplePDF.x);
+ −
+ − % Theor PDF
+ − theorPDF = copy(samplePDF,1);
+ − theorPDF.setY(Fpdf(samplePDF.x,dofX,dofY));
+ − theorPDF.setName('theor. PDF');
+ −
+ − % Theor PDF moments
+ − % theorK = utils.math.Kurt(samplePDF.x);
+ − % theorS = utils.math.Skew(samplePDF.x);
+ −
+ − % Plot PDFs
+ − if showPlots
+ − iplot(samplePDF,theorPDF);
+ − set(get(gca,'YLabel'),'String','PDF')
+ − set(get(gca,'XLabel'),'String','F statistic')
+ − end
+ −
+ − % Perform test on the hypothesis that both distributions should be
+ − % chi2-distributed or, equivalently, that the ratio of the two spectra
+ − % should be F-distributed. The comparison is actually done in chi2 sense.
+ − sampleHIST = hist(F,plist('N',nBin));
+ − ix = sampleHIST.y>=5;
+ − C = sum(sampleHIST.y)*mean(diff(sampleHIST.x));
+ − NN = sampleHIST.y(ix);
+ − nn = C*theorPDF.y(ix);
+ − r = (NN-nn).^2./nn;
+ − chi2 = sum(r);
+ − dof = numel(NN) - 1;
+ − Y2 = dof + sqrt(2*dof/(2*dof+sum(1./nn)))*(chi2-dof);
+ − critValue(1) = utils.math.Chi2inv(alpha/2,dof);
+ − critValue(2) = utils.math.Chi2inv(1-alpha/2,dof);
+ −
+ − % Perform the test on PDFs
+ − Chi2test = Y2<critValue(1) | Y2>critValue(2);
+ −
+ − % % Sample CDF
+ − % samplePDF = hist(F,plist('N',1000));
+ − % sampleCDF = copy(samplePDF,1);
+ − % sampleCDF.setY(cumsum(samplePDF.y)/sum(samplePDF.y));
+ − % sampleCDF.setName('sample CDF');
+ − %
+ − % % Theor CDF
+ − % theorCDF = copy(samplePDF,1);
+ − % theorCDF.setY(Fcdf(samplePDF.x,dofX,dofY));
+ − % theorCDF.setName('theoretical CDF');
+ − %
+ − % % Plot CDFs
+ − % iplot(sampleCDF,theorCDF);
+ − % set(get(gca,'YLabel'),'String','CDF')
+ − % set(get(gca,'XLabel'),'String','F statistic')
+ −
+ − % Critical values
+ − [test,critValue,pValue] = utils.math.Ftest(F.y,dofX,dofY,alpha,twoTailed);
+ −
+ − % Build AOs for critical values
+ − % if twoTailed
+ − critValueLB = ao(plist('xvals',X.x,'yvals',repmat(critValue(1),size(X.x)),...
+ − 'type','fsdata','fs',X.fs,'xunits',X.xunits,'name','critical values'));
+ − critValueUB = ao(plist('xvals',X.x,'yvals',repmat(critValue(2),size(X.x)),...
+ − 'type','fsdata','fs',X.fs,'xunits',X.xunits,'name','critical values'));
+ − % else
+ − % critValue = ao(plist('xvals',X.x,'yvals',repmat(critValue,size(X.x)),...
+ − % 'type','fsdata','fs',X.fs,'xunits',X.xunits,'name','critical values'));
+ − % end
+ −
+ − % % Build AOs for p-values and confidence levels
+ − % pValue = ao(plist('xvals',X.x,'yvals',pValue,...
+ − % 'type','fsdata','fs',X.fs,'xunits',X.xunits,'name','p-values'));
+ − % % pValue.setDy(pValueDy);
+ − %
+ − % if twoTailed
+ − % confLevelUB = ao(plist('xvals',X.x,'yvals',repmat(1-alpha/2/numel(F.y),size(X.x)),...
+ − % 'type','fsdata','fs',X.fs,'xunits',X.xunits,'name','confidence level'));
+ − % confLevelLB = ao(plist('xvals',X.x,'yvals',repmat(alpha/2/numel(F.y),size(X.x)),...
+ − % 'type','fsdata','fs',X.fs,'xunits',X.xunits,'name','confidence level'));
+ − % else
+ − % confLevel = ao(plist('xvals',X.x,'yvals',repmat(1-alpha/n,size(X.x)),...
+ − % 'type','fsdata','fs',X.fs,'xunits',X.xunits,'name','confidence level'));
+ − % end
+ −
+ − % Rejection index
+ − ix = find(test);
+ − if ~isempty(ix)
+ − H1 = F.setXY(plist('x',F.x(ix),'y',F.y(ix)));
+ − H1.setDy([]);
+ − H1.setName('H0 rejected');
+ − end
+ −
+ − % Compute the number of sigmal corresponding to the confidence level
+ − % Ns = erfinv(1-alpha)*sqrt(2);
+ −
+ − % Make the test: H0 rejected?
+ − % if twoTailed
+ − % % test = any( (F.y-Ns.*F.dy)>critValueUB.y | (F.y+Ns.*F.dy)<critValueLB.y );
+ − % test = any( F.y>critValueUB.y | F.y<critValueLB.y );
+ − % else
+ − % % test = any( (F.y-Ns.*F.dy)>critValue.y );
+ − % test = any( F.y>critValue.y );
+ − % end
+ −
+ − % Perform the F-test on spectra
+ − Ftest = any(test);
+ −
+ − % Plot results
+ − pl = plist('yscales',{'All', 'log'},'autoerrors',0);
+ − % if twoTailed
+ − if showPlots
+ − if Ftest
+ − [hfig, hax, hli] = iplot(F,critValueLB,critValueUB,H1,pl);
+ − set(hli(4),'linestyle','none','marker','s','MarkerSize',10,'MarkerEdgeColor','r');
+ − else
+ − [hfig, hax, hli] = iplot(F,critValueLB,critValueUB,pl);
+ − end
+ − set(hli(2:3),'color','r');
+ − % h2 = iplot(pValue,confLevelLB,confLevelUB,pl);
+ − % else
+ − % if Ftest
+ − % [hfig, hax, hli] = iplot(F,critValue,H1,pl);
+ − % set(hli(3),'linestyle','none','marker','s','MarkerSize',10,'MarkerEdgeColor','r');
+ − % else
+ − % [hfig, hax, hli] = iplot(F,critValue,pl);
+ − % end
+ − % set(hli(2),'color','r');
+ − % % h2 = iplot(pValue,confLevel,pl);
+ − % end
+ − set(get(gca,'YLabel'),'String','F statistic')
+ − end
+ −
+ − % Output final result
+ − test = any(Chi2test | Ftest);
+ −
+ − end
+ −