function [A,Sig]=FA(X,K,plotit,W)
format compact;
[D,N] = size(X);

if nargin < 3
   
   plotit=0;
     
 end

 [D,N] = size(X);
 Nshow = 100;
 threshold = 1e-3;
 Convergence=threshold+1;
 loop=0;
 Nloop=200;
 
   
 figure(31); clf;
 
% initialize
  Sig = eye(D);
  A = rand(D,K);
  
 % main loop
 while loop < Nloop & Convergence>threshold;
loop=loop+1;
    
   % E-step   
   IN = inv(A*A'+Sig);
   M = A'*IN * X;
   C = eye(K) - A'*IN*A;
   SC = N*C + M*M';
   
   % M-step   
   TA = (X*M') * inv(SC);
   TSig = 1e-5 + diag(diag(X*X' - TA*M*X'))/N;
   
   % Log-Likelihood   
   TL=0.5*(log(det(IN)) - sum(sum((IN*X).*X))/N );
   
   figure(31);
   plot(loop,TL,'go','Markersize',2,'Linewidth',2);hold on;
   drawnow;
   
   % Convergence
   if loop ~=1
      Convergence = TL - L;
      if Convergence < 0
         keyboard
      end
      
   end
  
   % Updates   
   A = TA;
   Sig = TSig;
   L=TL;
   
   %plotting stuff
if plotit == 1
      
if K==2
   
Axx=A*[-4;0];
Axy=A*[4;0];
Ayx=A*[0;-4];
Ayy=A*[0;4];
 
XX=[Axx(1) Axy(1)];
XY=[Axx(2) Axy(2)];
XZ=[Axx(3) Axy(3)];
YX=[Ayx(1) Ayy(1)];
YY=[Ayx(2) Ayy(2)];
YZ=[Ayx(3) Ayy(3)];

else
   
   Ax=A*-4;
   Ay=A*4;
   
   Xx=[Ax(1) Ay(1)];
   Yy=[Ax(2) Ay(2)];
   Zz=[Ax(3) Ay(3)];
end

figure(2);clf;
rotate3d on;
   plot3(X(1,:),X(2,:),X(3,:),'c+','Markersize',2,'Linewidth',2);hold on;
   if K==2
   plot3(W(1,:),W(2,:),W(3,:),'m','LineWidth',2);
   plot3(W(4,:),W(5,:),W(6,:),'m','LineWidth',2);
   plot3(XX,XY,XZ,'r','LineWidth',2);
   plot3(YX,YY,YZ,'r','LineWidth',2);
else
   plot3(W(1,:),W(2,:),W(3,:),'m','LineWidth',2);
   plot3(Xx,Yy,Zz,'r','LineWidth',2);
end
    drawnow;
end
   pause
end%while