% 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