Mercurial > hg > ltpda
view m-toolbox/classes/+utils/@math/SFtest.m @ 19:69e3d49b4b0c database-connection-manager
Update ltpda_uo.fromRepository
author | Daniele Nicolodi <nicolodi@science.unitn.it> |
---|---|
date | Mon, 05 Dec 2011 16:20:06 +0100 |
parents | f0afece42f48 |
children |
line wrap: on
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