function O = lyngby_fir_test_chi2(Y, X, B, varargin)

% lyngby_fir_test_chi2 - Chi square test for FIR filter
%
%       lyngby_fir_test_chi2(Y, X, B, 'PropertyName', 'PropertyValue')  
% 
%       Input:    Y   Data matrix 
%                 X   Design matrix
%                 B   Estimated filter
%
%       Property: RegMatrix  Regularization matrix
%
%       Output:   O   Lower tail P-values
%
%       A mass-univariate chi squared test with a null hypothesis of
%       B=0 for the FIR filter with the model: Y = X*B+U. The test is
%       performed separately for each column in B, ie, for each voxel
%       if Y is a (brain scans x voxels) data matrix.
%
%       The noise variance in U is estimated with the maximum
%       likelihood regarded as independent Gaussian distributed. When
%       the number of samples (eg, scans) is small and the noise
%       variance is estimated this chi2 test is biased, ie, is the
%       test is a large sample test.  
%
%       It is also possible to use this function to test a regularized
%       model, such as the 'smooth FIR' by inputting the
%       regulatization matrix with the 'RegMethod' property. This
%       introduces an other bias causing a too high estimate of the
%       variance. The two biases might cancel each other somewhat. 
%
%       Ref: Goutte, Nielsen, Hansen (2000), IEEE
%              Trans. Med. Imaging. 19(12):1188-1201, Equation 13. 
%
%       Example: 
%         XN = randn(70, 10000);     % eg, 70 scans and 10000 voxels
%         PN = randn(70, 1);
%         B = lyngby_fir_main(PN, XN);
%         G = lyngby_fir_paradigm2design(PN, size(B,1));
%         LT = lyngby_fir_test_chi2(XN, G, B);
%         hist(LT, 100), title('Histogram of lower tail area')
%
%       See also: LYNGBY, LYNGBY_FIR_MAIN,
%                 LYNGBY_FIR_PARADIGM2DESIGN, LYNGBY_FIR_REGMATRIX.
%
% $Id: lyngby_fir_test_chi2.m,v 1.4 2003/02/21 14:26:01 fnielsen Exp $


    % Default values
    sigma2 = [];
    dof    = [];
    R      = [];

    
    % Read properties
    n = 1;
    while n < length(varargin)
      arg = lower(varargin{n}); 

      if strcmp(arg, 'regmatrix')
	n = n + 1;
	arg = varargin{n};
	if isnumeric(arg) & size(arg,1) == size(arg,2) & ...
	      size(arg,1) == size(B,1)
	  R = arg;
	else
	  error(['The argument with ''RegMatrix'' should be a '...
	      'symmetric matrix']);
	end
	
      else
	error('Invalid property');
      end
      n = n + 1;
    end
    

    [N,P] = size(Y);
    Yhat = X * B;


    % Degrees of freedom
    if isempty(dof)
      dof = size(B,1);
    end
    
    % Noise mass-univariate variance sigma2 as a vector over P variables
    if isempty(sigma2)
      Uhat = Y - Yhat;
      sigma2 = sum(Uhat.^2,1) / N;
    end

    if isempty(R)
      % Without regularization
      chi2 = sum(Yhat.^2,1) ./ sigma2;
    else
      % With regularization
      for p = 1:P
	IW = R / sigma2(p);
	chi2(p) = B(:,p)' * (IW + X'*X/sigma2(p)) * B(:,p) / N;
      end
    end
    
    pvalue = lyngby_cdf_chi2(chi2, dof);
    
    O = pvalue;










