﻿// P81aMixIterSim.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);

nd=550;
// Simulace ----------------------------
thS=list();
thS(1)=[1 1]';
thS(2)=[2 5]';
thS(3)=[5 0]';
thS(4)=[8 3]';
thS(5)=[6 6]';
thS(6)=[4 3]';
nc=length(thS);
sdS=.6; cvS=[sdS**2 0;0 sdS**2];
alS=fnorm(1+ones(1,nc));
for t=1:nd
  cS(t)=sum(rand(1,1,'u')>cumsum(alS))+1;
  yt(:,t)=thS(cS(t))+sdS*rand(2,1,'n');
end // ---------------------------------

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

// Iterativní odhad
tht=[]; ct=[];
for x=1:50     // iterace
  for t=1:nd   // odhad pro všechna data
    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);                         // neuvažuje se alfa

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

  // chci, aby v komponentách byly maximální pravděpodobnosti
  // tedy ne: .3 .3 .4  ale: 0 0 1
  kr=0;
  for i=1:nc
    kr=kr+sum(wt(i,j(i).c));            // kriterium
  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                                   // pro stat. komponenty je střední
  clear wt;                             // hodnota průměr 
end

// Výsledky
[q T]=c2c(cS,c);                        // přetočení komponent
wrong=sum(q(c)(:)~=cS(:))               // počet špatných proti simulaci
from=nd

set(scf(1),'position',[500 100 600 500]);    // data a centra komponent
plot(yt(1,:),yt(2,:),'.','markersize',3)
for i=1:nc
  plot(thE(i)(1),thE(i)(2),'.','markersize',8)
end
