function GMM = GMMstat(data)

N = size(data,1);
u = mean(data);
w = (cov(data,1)+cov(data,1)')/2;
e = (data-u)*w^(-1/2);
k = size(data,2);

score = e.^3;
v_theta = 6*ones(k,1);

choose2 = nchoosek(1:k,2);
for i = 1:size(choose2,1)
    score = [score (e(:,choose2(i,1)).^2-1).*e(:,choose2(i,2))  (e(:,choose2(i,2)).^2-1).*e(:,choose2(i,1))];
end
v_theta = [v_theta; 2*ones(size(choose2,1)*2,1)];

choose3 = nchoosek(1:k,3);
for i = 1:size(choose3,1)
    score = [score e(:,choose3(i,1)).*e(:,choose3(i,2)).*e(:,choose3(i,3))];
end
v_theta = [v_theta; ones(size(choose3,1),1)];
v_theta = diag(v_theta);

GMM = N*mean(score)*v_theta^(-1)*mean(score)';

end

