// P71Mix1Stat.sce
// Mixture estimation - static components and pointer model
// - simulated data two-dimensional data
// - initialization by parameters from simulation + noise
// - with or withou estimation of covariances of components
//   FOR NC>5 IT IS NECESSARY TO SET PARAMETERS FOR SIMULATION
[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)

nd=100;             // number of data
nc=3;               // number of componentd
I_estCov=0;         // estimation of noise covariances ? 0|1 no|yes

// simulated reg.coef.
Sim.nc=nc;
Sim.Cy(1).th=[0.9 -.5]';
Sim.Cy(2).th=[0.4 0.6]';
Sim.Cy(3).th=[-.7 -.5]';
Sim.Cy(4).th=[-.1  .8]';
Sim.Cy(5).th=[ .5 -.1]';
// simulated noise covariances
Sim.Cy(1).sd=0.1*[1 0;0,1];
Sim.Cy(2).sd=0.1*[1 0;0,1];
Sim.Cy(3).sd=0.1*[1 0;0,1];
Sim.Cy(4).sd=0.1*[1 0;0,1];
Sim.Cy(5).sd=0.1*[1 0;0,1];
// simulated pointer parameters
Sim.Cp.th=fnorm(ones(1,nc)+.1,2);

Sim.ct(1)=1;                  // initial poiner

// initial parameters
a=.8;                         // std of scattering initial parametrs
for j=1:nc                    //  from those used in simulation
  [mr,mc]=size(Sim.Cy(j).th);
  Ps=[Sim.Cy(j).th;1]+[a*rand(mr,mc,'n');0]; // initial parameters
  Est.Cy(j).V=Ps*Ps';                        // statistics
  Est.Cy(j).th=Sim.Cy(j).th+a*rand(2,1,'n'); // pt.est. of reg.coef
  Est.Cy(j).sd=.1*eye(2,2);                  // standard deviation
end
Est.ka=ones(1,nc);            // counter
Est.Cp.V=ones(1,nc);          // pointer statistics
Est.Cp.th=fnorm(ones(1,nc));  // pointer parameter
w=fnorm(ones(1,nc));          // weights

// SIMULATION ==========================================================
for t=2:nd
  Sim.ct(t)=sum(rand(1,1,'u')>cumsum(Sim.Cp.th))+1;    // pointer
  j=Sim.ct(t); // active component
  y=Sim.Cy(j).th+Sim.Cy(j).sd*rand(2,1,'norm');        // output
  Sim.yt(:,t)=y;
end   

// ESTIMATION ==========================================================
printf(' ')
for t=2:nd  
  if t/10==fix(t/10), printf('.'); end
  for j=1:nc
    [xxx,G(j)]=GaussN(Sim.yt(:,t),Est.Cy(j).th,Est.Cy(j).sd); 
                                             // proximity
  end
  Lq=G-max(G);
  q=exp(Lq);

  ww=q'.*Est.Cp.th; w=ww/sum(ww);            // generation of weights
  wt(:,t)=w';

  // Update of statistic 
  Ps=[Sim.yt(:,t)' 1];                       // extended reg.vec.
  for i=1:nc
    Est.Cy(i).V=Est.Cy(i).V+w(i)*Ps'*Ps;     // information matrix
    Est.ka(i)=Est.ka(i)+w(i);                // counter
    Est.Cp.V(i)=Est.Cp.V(i)+w(i);            // pointer statistics

    //nove rozdeleni informacni matice V
    Vyy=Est.Cy(i).V(1:2,1:2);                // part Vyy - psi.psi'
    Vy=Est.Cy(i).V($,1:2);                   // part Vy  - psi.y
    V1=Est.Cy(i).V($,$);                     // part V1  - y.y
    Est.Cy(i).th=inv(V1+1e-8*eye(V1))*Vy;    // pt.est. - reg.coef.
    Est.Cy(i).tht(:,t)=Est.Cy(i).th';        // pt.est. - covar.
    if I_estCov~=0
      // pt.est. of noise covariance - used or not
      Est.Cy(i).cv=(Vyy-Vy'*inv(V1+1e-8*eye(V1))*Vy)/Est.ka(i); 
    end
  end
  Est.Cp.th=fnorm(Est.Cp.V,2);     // pt.est. of pointer parameter
  [ss,Est.ct(1,t)]=max(w);         // store
end

// Results
disp(Sim.Cp.th,'pt.pars_sim')
disp(Est.Cp.th,'pt.pars_est')

s=2:nd;
wr=sum(Sim.ct(s)~=Est.ct(s)');
printf('\n Wrong classifications %d from %d\n',wr,length(s))

set(scf(1),'position',[50 50 400 400])
plot(s,Sim.ct(s),'bo',s,Est.ct(s),'rx')
legend('simulated','estimated');
title 'Values of the pointer'

set(scf(2),'position',[550 50 400 600])
for i=1:nc
  subplot(nc,1,i)
  plot(Est.Cy(i).tht')
  title('Parameters estimate evolution '+string(i))
end

