function [StatResuV,  StatResuP, StatResuPc, MVResu] = mars_stat_compute(...
				xCon, Xs, V, betas, RMS);
%--------------------------------------------------------------------
%
%--------------------------------------------------------------------
nROI = size(betas,2);
nCon = length(xCon);

xpxm = spm_sp('xpx-',Xs);
xpVx = Xs.X'*V*Xs.X;
Bcov = xpxm*xpVx*xpxm;

[trRV trRVRV] = spm_SpUtil('trRV',Xs,V);
erdf = trRV^2/trRVRV;


T_indices = [];
F_indices = [];
dfnum     = [];

StatResuV    = zeros(nCon, nROI);
StatResuP    = zeros(nCon, nROI);
check_Tvalue = zeros(nCon, nROI);
MVResu       = zeros(nCon, 1);

for ii = 1:nCon

%	[edf_tsp edf_Xsp] = spm_FcUtil('FconEdf', xCon(ii), Xs, V);

	switch(xCon(ii).STAT)
	   case 'T'
		%----------- to check calculation with h -----------
		h       = spm_FcUtil('Hsqr', xCon(ii), Xs);
		[trMV trMVMV]= spm_SpUtil('trMV', ...
				spm_FcUtil('X1o',xCon(ii),Xs), V);

		check_Tvalue(ii,:) = ((h*betas)/trMV)./RMS;
		check_dfnum(ii) = trMV;

		%----------- t value  -----------

    		StatResuV(ii,:) = (xCon(ii).c'*betas) ./ ...
                                   (RMS .* sqrt((xCon(ii).c'*Bcov*xCon(ii).c)));
		StatResuP(ii,:) = 1 - spm_Tcdf(StatResuV(ii,:), erdf);
		T_indices = [T_indices ii];


    	   case 'F'  
		[trMV trMVMV]= spm_SpUtil('trMV', ...
				spm_FcUtil('X1o',xCon(ii),Xs), V);
		dfnum   = [dfnum trMV^2/trMVMV];
		h       = spm_FcUtil('Hsqr', xCon(ii), Xs);

		StatResuV(ii,:) = (sum( (h*betas).^2 )/trMV) ./ (RMS.^2);
		check_Tvalue(ii,:) = StatResuV(ii,:) ;
		StatResuP(ii,:) = 1 - spm_Fcdf(StatResuV(ii,:), ...
					[dfnum(end) erdf]);
		F_indices       = [F_indices ii];
		

	   otherwise
	        error(['unknown STAT "',xCon(ii).STAT,'"'])
	end
end

%- corrected P value for the number of ROI
%--------------------------------------------------------------------	   
StatResuPc = ones(nCon, nROI) - (ones(nCon, nROI) - StatResuP).^nROI;



