﻿// P75aMixExp.sce
// Pokusy s exp proximity jako likelihood
// - 3 komponenty
// - 2D data
// ------------------------------------------
[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

nix=10;
Sim.Cx(1).th=[.5 5]'; 
Sim.Cx(2).th=[6 .5]';
Sim.Cx(3).th=[13 15]';
Est.Cpx.th=[.5 .2 .3];
Est.Cpx.V=.1*ones(1,3);                     // pointer statistics
Est.kax=nix*ones(1,3);                             // counter
cx(1)=1;

nd=1000;
for t=2:nd
  cx(t)=sum(rand(1,1,'u')>cumsum(Est.Cpx.th))+1;
  u=rand(2,1,'u'); a=Sim.Cx(cx(t)).th;
  x(:,t)=-log(u).*a;
end

Vx=[Sim.Cx(1).th Sim.Cx(2).th Sim.Cx(3).th]; 
for i=1:3
  Est.Cx(i).V=Vx(:,i)*nix;
  Est.Cx(i).th=Est.Cx(i).V;
end

printf(' ')  //  ============= TIME LOOP ==============
for t=2:nd  
  if t/fix(nd/10)==fix(t/fix(nd/10)), printf('.'); end
  
  // váhy 1. ÚSEK
  for j=1:3
    [xxx,Gx(j)]=ExpN(x(:,t),Est.Cx(j).th); // model
//    [xxx,Gx(j)]=ExpMLik(x(:,t),Est.Cx(j).V,Est.kax(j));  // likelihood
//    [xxx,Gx(j)]=GaussN(x(:,t),Est.Cx(j).th,.01*eye(2,2));
//      dx=x(:,t)-Est.Cx(j).th; dx=dx(:);
//      Gx(j)=1/(dx'*dx)**5;
  end

  Lx=Gx-max(Gx);              // rough normlization
  qx=exp(Lx);                 // exponent

  ww=qx'.*Est.Cpx.th;         // matrix weights
  wx=ww/sum(ww);              // normalization    - f(c(t),c(t-1)|d(t))
  wtx(:,t)=wx';               // stor
//xx=x(:,t)', wx,cxx=cx(t)
//if t>10, pause, end
  // Update of statistic 1ST COMPONENT ==================
  Est.kax=Est.kax+wx;                           // counter
  Est.Cpx.V=Est.Cpx.V+wx;                       // ptr.stat. update
  Psx=[x(:,t)' 1];                       // extended reg.vec.
  for i=1:3
    // přepočet statistiky
    Est.Cx(i).V = Est.Cx(i).V + wx(i)*x(:,t);     // information matrix

    //bodové odhady
    Est.Cx(i).th=Est.Cx(i).V./(ones(2,1)*Est.kax(i));
    Est.Cpx.th=fnorm(Est.Cpx.V);
  end
end
disp ' '
[xxx,Ecx]=max(wtx,'r');

  //USSAGE:       ct=q(Ect)
  // set:  [q T]=c2c(ct,Ect);      simul and estim pointer
  // plot(1:nd,ct,1:nd,q(Ect))     plotting 

[q,T]=c2c(cx',Ecx);
Ecc=q(Ecx);
wrong=sum(cx'~=Ecc')
from=nd

nc=3;
// sim -----------------------
CS=list();
for i=1:nc
  j=find(cx==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
  if ~isempty(C(i))
  plot(C(i)(1,:),C(i)(2,:),tx(i),'markersize',3)
  end
end
title 'Estimated clusters'
set(gca(),'data_bounds',[0 45 0 40])
