﻿// P75bMixGamma.sce
// Pokusy s gamm proximity jako likelihood
// - dvourozměrná data
// - více komponent
// - dvourozměrna 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
rand('seed',0);

nc=5;          // number of components
nv=2;          // number of variables
nd=1000;       // number of data
prox=1;        // 1=gamm proximity, 0=gaussian-like prox.

thS=list();
//      a b
thS(1)=[1 2      // x1
       10 2];    // x2        Egam=a/b 
thS(2)=[8 1
        1 2];
thS(3)=[8 2
       12 3];
thS(4)=[28 2
       5 3];
thS(5)=[.8 2
        .2 3];
alS=[.2 .1 .3 .2 .2]; alS=fnorm(alS(1:nc));
cS(1)=1;

for t=2:nd     // ============== SIMULATION ============
  cS(1,t)=sum(rand(1,1,'u')>cumsum(alS))+1;
  u=rand(nv,1,'u'); 
  thP=thS(cS(t));
  x(:,t)=cdfgam("X",thP(:,1),thP(:,2),u(:),1-u(:));
end
//scatt(x)
//return

              // ============ INITIALIZATION ============
ni=50;
ka=ni*(.1*rand(1,nc,'u')+1);
nu=ni*(rand(1,nc,'u')+5);          // zatim se nepoužije (s tím je horší)
al=nu/sum(nu);                     // \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
zz=[];
for i=1:nc
  z=zeros(2,ni);
  for j=1:ni
    u=rand(nv,1,'u');
    P=thS(i); //+2*rand(thS(i),'u'); 
    P(:,2)=ones(2,1);
    z(:,j)=cdfgam("X",P(:,1),P(:,2),u,1-u);
  end
  Cy.V=[sum(z,2) sum(log(z),2)];
  zz=[zz z];

  s=log(Cy.V(:,1)/ka(i))-Cy.V(:,2)/ka(i);
  Cy.th(:,1)=(3-s+sqrt((3-s)^2+24*s))./(12*s);
  Cy.th(:,2)=Cy.th(:,1)*ka(i)./Cy.V(:,1);    

  Est.Cy(i)=Cy;
end
ct=ones(1,nd); 
tht=list();
for i=1:nc, tht(i)=zeros(4,nd); end

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

  // proximity
  for i=1:nc
    if prox==1
      [xxx,G(i)]=GammaG(x(:,t),Est.Cy(i).th); 
    else
      Ex=Est.Cy(i).th(:,1)./Est.Cy(i).th(:,2);
      G(i)=(x(:,t)-Ex)'*(x(:,t)-Ex)/1000;
    end    
  end
  Gn=G-max(G);
  g=exp(Gn);
//  g=g.*al';
  w=g/sum(g);
  wt(:,t)=w(:);

  // přepočet statistik pro gamma     
  for i=1:nc
    Est.Cy(i).V(:,1)=Est.Cy(i).V(:,1)+w(i)*x(:,t);        // data
    Est.Cy(i).V(:,2)=Est.Cy(i).V(:,2)+w(i)*log(x(:,t));   // log()   

    // bodové odhady parametrů komponent pro gamma
    s=log(Est.Cy(i).V(:,1)/ka(i))-Est.Cy(i).V(:,2)/ka(i);
    Est.Cy(i).th(:,1)=(3-s+sqrt((3-s)^2+24*s))./(12*s);
    Est.Cy(i).th(:,2)=(Est.Cy(i).th(:,1)*ka(i))./Est.Cy(i).V(:,1);   
    tht(i)(:,t)=Est.Cy(i).th(:);
  end
  ka=ka+w';
  nu=nu+w';    // zatim se nepoužije
  al=nu/sum(nu);
end
[xxx,cE]=max(wt,'r');
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.';'m.';'k.';'y.'];
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])

set(scf(3),'position',[10 50 1500 400])
for i=1:nc
subplot(1,nc,i),plot(tht(i)')
end
legend('11','21','12','22');

