function factorPCA(X,names);
%function factorPCA(X,names)
% 
% X = variabili
% names = etichette (omettere la seconda variabile names se non usata)
%

[n,p]=size(X);
%if (nargin<2),names=[];for k=1:p,names=[names;num2str(k)];end,end
if (nargin<2), names=int2str((1:p)'); end
m=mean(X);
X1=X-repmat(m,n,1);
R=corrcoef(X1);
fprintf('Matrice di correlazione\n')
disp(R)

[V,E]=eig(R);E=abs(real(E));V=real(V);

[e,id]=sort(diag(E));
e=e(p:-1:1);
V=V(:,id(p:-1:1));

fprintf('eigs e Varianza spiegata da ogni componente\n')
disp([e,e/sum(e)])
fprintf('Varianza cumulata \n')
disp(cumsum(e)./sum(e))

fprintf('Intervalli di confidenza 0.95 per la varianza spiegata \n')
disp([e./(1-norminv(0.05/2)*sqrt(2/n)), e./(1+norminv(0.05/2)*sqrt(2/n))]./p)

%figure(1)
subplot(2,2,3);
plot(e,'*-');hold on;plot(cumsum(e)/sum(e),'r*-');hold off
title('scree plot');xlabel('# componente principale')
ylabel('varianza spiegata')


L1=V*sqrt(diag(e));
fprintf('carichi \n')
disp(L1)

%k=1
fprintf('\n')
fprintf('1 carico\n')
S1=L1(:,1)*L1(:,1)';
psi=diag(diag(R-S1));
fprintf('Errore nella matrice di corr.\n')
disp(R-(S1+psi))
fprintf('Norma-quadro dell''errore nella matrice di corr.\n')
disp(norm(R-(S1+psi), 'fro')^2)

%k=2
fprintf('\n')
fprintf('2 carichi\n')
S1=L1(:,1:2)*L1(:,1:2)';
psi=diag(diag(R-S1));
fprintf('Errore nella matrice di corr.\n')
disp(R-(S1+psi))
fprintf('Norma-quadro dell''errore nella matrice di corr.\n')
disp(norm(R-(S1+psi), 'fro')^2)

%figure(2)
subplot(2,2,1);
biplot(L1(:,1:2),'VarLabels',names)
title('carichi senza rotazione')

%figure(3);
subplot(2,2,2);
[Lr,TT]=rotatefactors(L1(:,1:2));
biplot(Lr(:,1:2),'VarLabels',names)
title('carichi dopo rotazione varimax')

fprintf('\n')
fprintf('2 Carichi ruotati\n')
disp(Lr(:,1:2))

s=std(X); X1=(X-ones(n,1)*mean(X))/diag(s);

fprintf('\n')
fprintf('Stima fattori \n')
osserv=1;
(Lr'*Lr)\(Lr'*(X1(osserv,:)-mean(X1))')

%k=3
fprintf('\n')
%figure(4);
subplot(2,2,4);
Lr3=rotatefactors(L1(:,1:3));
fprintf('3 carichi ruotati\n')
disp(Lr3)
biplot(Lr3,'VarLabels',names)
S1=Lr3*Lr3';
psi=diag(diag(R-S1));
fprintf('Errore nella matrice di corr.\n')
disp(R-(S1+psi))
fprintf('Norma-quadro dell''errore nella matrice di corr.\n')
disp(norm(R-(S1+psi), 'fro')^2)

fprintf('\n')
fprintf('Stima fattori \n')
(Lr3'*Lr3)\(Lr3'*(X1(1,:)-mean(X1))')
(L1(:,1:3)'*L1(:,1:3))\(L1(:,1:3)'*(X1(1,:)-mean(X1))')

