```% Calculates a measure of multivariate association of Gaier, 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:
% Gaier, S., Ruppert, M., & Schmid, F. (2010). A multivariate version of Hoeffdings Phi-Square. Journal of Multivariate Analysis, 101(10), 25712586.

% 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

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```