
% Calculates a measure of multivariate association of Gaißer, S., Ruppert, M., & Schmid, F. (2010)
% Inputs: data - n x d matrix containing n realizations of d random
% variables, association between which is to be measures.
% Output: phi - 1x1 measure of association (phi = 0 corresponds to mutual
% independence of variables in columns of data)

% References:
% Gaißer, S., Ruppert, M., & Schmid, F. (2010). A multivariate version of Hoeffding’s Phi-Square. Journal of Multivariate Analysis, 101(10), 2571–2586.

% Notice: Code (C) 2014 Ivan Medovikov. Free for non-commercial use. No warranty, express or implied.
% Code may not be distributed without this notice. See http://medovikov.me
% for latest version. Contact: imedovikov@brocku.ca


function phi = multphi2(data)

n = size(data,1);
d = size(data,2);

%% Calculate empirical marginal d.f.-s

ranks = zeros(n,d);

for i=1:d
    for j=1:n
        ranks(j,i) = (1/n)*sum(data(:,i)<=data(j,i));
    end
end

%% Calculate multivariate Phi2

first = 0;
second = 0;

% First term of the sum

for j=1:n
    for k=1:n
        inner_product = 1;
        for i=1:d
            inner_product = inner_product * (1-max(ranks(j,i),ranks(k,i)));
        end
        first = first + inner_product;
    end
end

% Second term of the sum

for j=1:n
    inner_product = 1;
    for i=1:d
        inner_product = inner_product * (1-ranks(j,i)^2);
    end
    second = second + inner_product;
end

% Final statistic
       
numerator = ((1/n)^2) * first - ((2/n)*(1/2)^d) * second + (1/3)^d;

denominator_product = 1;
for i=0:d
    denominator_product = denominator_product * (i + 0.5);
end

denominator = (2/((d+1)*(d+2)))...
    - (1/(2^d))*(factorial(d)/denominator_product) + (1/3)^d;

phi = numerator / denominator;

end