// P71Mix3Pred.sce
// Prediction with known mixture (without estimation)
// - static components
// - dynamic pointer model
// - simulated one-dimensional data
// - known model parameters
[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                 

np=5;                              // number of steps of prediction
nc=3;                              // number of components 
nd=np+1+1000;                      // number of data 

// INITIALIZATION of simulation --------------------------------------
thS=[-2 1 5];                      // centres of components
cvS=[.1 .1 .1]*1; 
select 1       // these pointr models are pre-set, you can either
               // chose one or write your own model. Rows must sum to 1
case 1 then, alS=[0 1 0; 0 0 1; 1 0 0];
case 2 then, alS=[.07 .87 .06; .06 .07 .87; .92 .04 .04];
case 3 then, alS=[.8 .1 .1;.1 .2 .7;.1 .1 .8];  // parameters of the pointer 
end 
c(1)=1;                            // initial values 

// SIMULATION
for t=2:nd
  rnd=rand(1,1,'unif');
  cus=cumsum(alS(c(t-1),:));
  c(t)=sum(rnd>cus)+1;             // pointer generation 
  th=thS(c(t));
  sd=sqrt(cvS(c(t)));
  y(t)=th+sd*rand(1,1,'norm');     // output generation 
end

w1=ones(nc,1)/nc;                  // initial pointer distribution 

// PREDICTION
for t=2:nd // TIME LOOP ----------------------------------------------
  yt=y(t);                         // measured output

  // proximity of yt to components
  for i=1:nc
    m(i)=GaussN(yt,thS(i),cvS(i)); // density of normal component
  end  
  
  // weighting vector 
  wP=fnorm(alS'*w1);               // weights for prediction
  wp=(alS')^np*wP;                 // np-steps prediction
  w=fnorm(wP.*m);                  // normalization
  
  // update of statistics
  Ps=[yt;1];                       // ext. regression vector 
  yp=0;
  for i=1:nc
    yp=yp+wp(i)*(thS(i)+sqrt(cvS(i))*rand(1,1,'norm'));
  end
  ytp(t+np)=yp;
 
  w1=w;                            // old weighting vector 
  wt(:,t)=w;                       // stor weighting vector 
  wtp(:,t+np)=wp;
end  // END OF TIME LOOP ---------------------------------------------

// RESULTS
[xxx ce]=max(wtp,'r');              // estimated active components
s=(np+2):nd;
wr=sum(c(s)~=ce(s)'); 
tx='\nWrong class. %d from %d, i.e. %d percent\n';
printf(tx,wr,length(s),100*wr/length(s))
s=(nd-100+1):nd;
plot(s,c(s),'bo','markersize',8)
plot(s,ce(s),'r.','markersize',3)
title 'Pointer and its estimate'
legend('pointer','estimate');
set(gca(),'data_bounds',[min(s) max(s) .9 nc+.1])
set(gcf(),'figure_position',[400,100])

scf(1);
set(gcf(),'figure_position',[600,150])
title 'Histograms of data and predicted data'
subplot(211)
[f1,n1]=histplot(150,y);
ylabel 'data'
subplot(212)
[f2,n2]=histplot(150,ytp);
ylabel 'predictions'

s=551:600;
scf(2);
set(gcf(),'figure_position',[200,150])
title 'time course of data and predictions'
plot(s,y(s),'ob')
plot(s,ytp(s),'xr')
legend('data','predictions');
