﻿// P81bMixIterRe.sce
// Pokus o odhad směsi iterativní off-line metodou
// - pro dané hodnoty parametrů se spočtou všechny váhy
// - z vah se určí celé ukazovátko 
// - pro něj se udělá nový odhad parametrů
// -- data se rozdělí podle komponent
// -- pro každou komponentu se spočte střední hiodnota
// - vypočte se kriterium (součet vah aktivních komponent) - není třeba
// - a to se opakuje
exec SCIHOME/ScIntro.sce, mode(0), rand('seed',0);

load _data/d_con.dat d_con;
yt=d_con';
[nv,nd]=size(yt);

nc=3;

// Inicializace
thI=list(); thE=thI;
for i=1:nc
  thI(i)=rand(nv,1,'u');
end
thE=thI;
cv=.1*eye(nv,nv);

tht=[]; ct=[];
for x=1:50     // iterace
  for t=1:nd   // off-line odhad
    for i=1:nc                     // proximity
      [xxx,G(i)]=GaussN(yt(:,t),thE(i),cv);
    end  
    Lq=G-max(G);                   // rough normlization
    w=exp(Lq);                     // exponent   
    w=w/sum(w);                              // není alfa

    wt(:,t)=w;
  end
  [xxx,c]=max(wt,'r');             // nové c
  if prod(c==ct($,:)), break, end  // test na změnu céček
  for i=1:nc
    j(i).c=find(c==i);             // komponenty
  end
  ct=[ct; c];                      // stará + nová   céčka

  kr=0;
  for i=1:nc
    kr=kr+sum(wt(i,j(i).c));       // kriterium (váhy -> 0 nebo 1)
  end
  printf('iter = %d\n',x)                      // iter
  printf('krit = %g\n\n',kr)                   // krit

  for i=1:nc
    thE(i)=mean(yt(:,j(i).c),2);   // nový odhad komponent
  end
  clear wt;
end

// Výsledky
// Dále je možno volit i1 a i2 (veličiny) a spouštět jako blok
cla
i1=2; i2=5;
scf(1);
plot(yt(i1,:),yt(i2,:),'.','markersize',3)
for i=1:nc
  plot(thE(i)(i1),thE(i)(i2),'.r','markersize',8)
end
