function factorMV(X,m,names);
%function factorMV(X,m,names);

stats.chisq=[];
format short
[n,p]=size(X);
%if (nargin<3),names=[];for k=1:p,names=[names,int2str(k)];end,end
if (nargin<3),names=[];for k=1:p,names=[names;k];end, names=int2str(names); end

R=corrcoef(X);
S=cov(X);
[lambda, psi, t, stats]=factoran(X,m);
lambda1=lambda;
lambda=lambda*t';
fprintf('\n')
fprintf('Carichi non ruotati\n')
disp(lambda)
fprintf('\n')
fprintf('Carichi ruotati\n')
disp(lambda1)

if (m>1 & m <4)
subplot(1,2,1); biplot(lambda(:,1:m),'VarLabels',names)
title('carichi senza rotazione')
subplot(1,2,2); biplot(lambda1(:,1:m),'VarLabels',names)
title('carichi dopo la rotazione')
end


fprintf('\n')
fprintf('Stime delle Varianze specifiche\n')
disp(psi)

fprintf('\n')
fprintf('Matrice errore \n')
disp(R-(lambda*lambda'+diag(psi)))

fprintf('\n')
fprintf('Varianza spiegata \n')
fprintf('Factor n.   proportion   cum.proport. \n')
vari=(sum(lambda.*lambda,1)/p)';
disp([(1:m)',vari, cumsum(vari)])

fprintf('\n')
fprintf('Ipotesi nulla sul numero di carichi (alpha=0.05)\n')

if isfield(stats,'chisq'),
fprintf('\n')
disp([stats.chisq,chi2inv(.95,stats.dfe)])
fprintf('\n')
if (stats.chisq>chi2inv(0.95,stats.dfe)),
  fprintf('reject H0\n')
else
  fprintf('do NOT reject H0\n')
end
else
fprintf('no test di ipotesi automatico\n')

coef=det(lambda1*lambda1'+diag(psi))/det(S);
st_new=(n-1-(2*p+4*m+5)/6)*log(coef);
chi_ref= chi2inv(.95, ( (p-m)^2-p-m)/2);
%disp([st_new,chi_ref])
if (st_new<0), fprintf('chi2 negativo\n');end
if (st_new>chi_ref), 
fprintf('numero di carichi non adatto\n');elseif (st_new>0)
fprintf('numero di carichi adatto\n');
end
end
