% Script to calculate Hotelling's T^2 statistic for comparing
% the means of multivatiate data.
% ---------------------
% Used for comparing two clusters of foci, to verify that they 
% are different.  The input files with coordinates should be
% a matrix of (num_of_voxels X [x y z]).
% ============================================================
% Kalina Christoff                                  2000-02-06



% clustA=load('2rel_1rel.group_1.roi');

fprintf('\nPlease, enter the names of two text files, one for each\n');
fprintf('cluster. The files should contain the x,y,z coordinates for\n');
fprintf('each voxel. Rows should be voxels.\n\n');

clustA = input('Name of file for cluster A: ','s');
clustB = input('Name of file for cluster B: ','s');

clustA = load(clustA);
clustB = load(clustB);

% calculate number of coordinates in each cluster
N_clustA   = size(clustA,1);
N_clustB    = size(clustB,1);

% calculate the vector of sample means for each cluster
M_clustA   = mean(clustA);
M_clustB    = mean(clustB);

% calculate the variance-covariance matrices for each cluster
COV_clustA = cov(clustA);
COV_clustB  = cov(clustB);


% calculate the pooled estimate of the variance-convariance matrix
% (assuming the population covariance matrices are the same)
% COV = {(n1-1)COV1 + (n2-1)COV2} / (n1+n2-2)

COVr_pooled = ( (N_clustA - 1)*COV_clustA  +  (N_clustB - 1)*COV_clustB )...
                    / ( N_clustA + N_clustB - 2);


% compute Hotelling's T-square statistic
% T^2 = n1*n2* ([x_bar1] - [x_bar2]) * COV^(-1) * ([x_bar1] - [x_bar2])' /
%                          / (n1+n2)

hotellingT = N_clustA*N_clustB*(M_clustA-M_clustB)*inv(COVr_pooled)*...
                  (M_clustA-M_clustB)' / (N_clustA+N_clustB);

% Hotelling's T can be conversted into a statistic that follows the
% F-distribution, with df=[num_of_variables, (n1+n2 - num_of_variables -1):
% F = (n1+n2 - num_of_variables -1) * T^2 / {(n1+n2 -2)* num_of_variabls}

F  = (N_clustA+N_clustB -2 -1) * hotellingT / ( (N_clustA+N_clustB -2) * 2 )
df = [2, (N_clustA + N_clustB -2 -1)]
p  = fpdf(F,df(1),df(2))

