﻿// P75cMixRay.sce
// Pokusy s exp proximity jako likelihood
// - dvourozměrná data
// - 3 komponenty
// ------------------------------------------
[u,t,n]=file();                         // find working directory
chdir(dirname(n(2)));                   // set working directory
clear("u","t","n")                      // clear auxiliary data
exec("ScIntro.sce",-1),mode(0)          // intro to sesion

nc=3;
nx=2;
nd=500;
rS=[.5 100 40
     50 130 1];
alS=[.3 .4 .3];
cS(1)=1;

for t=2:nd     // ============== SIMULATION ============
  cS(1,t)=sum(rand(1,1,'u')>cumsum(alS))+1;
  u=rand(nx,1,'u'); 
  x(:,t)=sqrt(-2*rS(:,cS(t)).*log(u));
end
//scatt(x)
//return
V=10*rand(2,3,'u')*10;
ka=[1 1 1]*10;
nu=[.5 .7 .6];
r=[.5 15 55
     1 30 105];

printf(' ')  //  ============= TIME LOOP ==============
for t=2:nd  
  if t/fix(nd/10)==fix(t/fix(nd/10)), printf('.'); end

  for i=1:nc
    G(1,i)=1;
    for j=1:nx
      G(1,i)=G(1,i)*log(x(j,t))-log(r(j,i))-.5*x(j,t)**2./r(j,i);
    end
  end
  Gn=G-max(G);
  g=exp(Gn);
  w=g/sum(g);
  wt(:,t)=w(:);

  // Update of statistic
  for i=1:nc
    V(:,i)=V(:,i)+w(i)*x(:,t)**2;
  end
  ka=ka+w;
  nu=nu+w;
 
  // Bodové odhady
  for i=1:nc
    r(:,i)=.5*V(:,i)/ka(i);
  end

end
printf('\n')

// VÝSLEDKY
[xxx,Ect]=max(wt,'r');
[q,T]=c2c(cS,Ect);
Ecc=q(Ect);
wrong=sum(cS(:)~=Ecc(:))
from=nd
 
// sim -----------------------
CS=list();
for i=1:nc
  j=find(cS==i);
  CS(i)=x(:,j);
end

tx=['b.';'r.';'g.'];
scf(1);
for i=1:nc
  plot(CS(i)(1,:),CS(i)(2,:),tx(i),'markersize',3)
end
title 'Simulated clusters'
set(gca(),'data_bounds',[0 45 0 40])

// est ------------------------
C=list();
for i=1:nc
  j=find(Ecc==i);
  C(i)=x(:,j);
end

set(scf(2),'figure_position',[850 200]);
for i=1:nc
  plot(C(i)(1,:),C(i)(2,:),tx(i),'markersize',3)
end
title 'Estimated clusters'
set(gca(),'data_bounds',[0 45 0 40])
