﻿// P74MixUnif.sce
// Mixture estimation with uniform components
//  Data:  cons. speed. gas. moment revs. slope
// -------------------------------------------
[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

function unifplot(th,s)
  // plot of a rectangle 
  // th   [L1 U1
  //       L2 U2]
  // s    format of plot

  if argn(2)<2, s='b'; end
  
  for j=1:max(size(th))
    L1=th(j)(1,1);
    R1=th(j)(1,2);
    L2=th(j)(2,1);
    R2=th(j)(2,2);

    plot([L1,R1],[L2,L2],s)
    plot([L1,R1],[R2,R2],s)
    plot([L1,L1],[L2,R2],s)
    plot([R1,R1],[L2,R2],s)
  end
endfunction

// Real data: individual variables (vith denotation of their meaning)
// 1  SPOTREBA      y1         9  PRIC-ZRYCH    v3        y - output
// 2  SPEED (2+3)/2 y2        10  STAC-SPEED    v4        u - control
// 3  VOLANT        u1        11  ROAD-CAN      v5        v - ext. input
// 4  PLYN          u2        12  CAS-UTC       v6
// 5  BRZDA         u3        13  UTM-X2        v7
// 6  R-STUP-A      u4        14  UTM-Y2        v8
// 7  MOMENT-MOT    v1        15  ALTITUDE2     v9
// 8  OTACKY-MOT    v2        16  COURSE       v10
// -----------------------------------------------------------------------------

nd=6000;                      // number of data (max 75420/5 = 15084)

// DATA LOAD
load _data/dataKF.dat
dt=DataKF.dataOrig;           // data from driven car
dt=dt([2 7],1:5:$);           // each 5th sample is taken
                              // selected data are: speed and moment of engine
// INITIALIZATION
y=dt(:,1:nd);
th=list();                    // initial parameters (lower and upper bounds)
th(1)=[60 70
       50 100];
th(2)=[60 70
       0 5];
th(3)=[60 70
       -22 -18];
th(4)=[130 140
       100 110];
th(5)=[120 130
       0 5];
thI=th;                       // remember initial setting
ny=size(th(1),1);             // dimension of data (given by th)
nc=length(th);                // number of components (given by th)
nu=ones(1,nc);                // pointer statistics
al=fnorm(nu);                 // pointer parameter
m=zeros(1,nc);
Lm=zeros(1,nc);
E=list(); D=list();
thL=list(); thU=list();       // upper and lower borders
for j=1:nc
  thL(j)=zeros(ny,nd);
  thU(j)=zeros(ny,nd);
end  
if ny~=size(dt,1), disp('ERROR: Wrong number of variables'), return, end // check

printf(' '), tt=fix(nd/10); 
for t=(2:nd)   // ----------- TIME LOOP --------------
  if t/tt==fix(t/tt), printf('.'); end
  
  // weights
  for j=1:nc
    E(j)=(th(j)(:,2)+th(j)(:,1))/2;          // expectation
    D(j)=(th(j)(:,2)-th(j)(:,1))**2/12;      // variance
    [xxx,Lm(j)]=GaussN(y(:,t),E(j),diag(D(j)));
  end
  Lm=Lm-max(Lm);
  m=exp(Lm);
  if sum(m)<-1E10, m=ones(1,nc); end
  w=m/sum(m);
  wt(:,t)=w'; 
     
  // on-line estimation - statistics
  if 0, nn=ones(1,nc); else nn=nu; end  // forgetting: yes or not
  for j=1:nc
    for i=1:ny
      yL=y(i,t)-th(j)(i,1);
      if yL<=0, th(j)(i,1)=th(j)(i,1)+w(j)*yL/nn(j); end
      yU=y(i,t)-th(j)(i,2);
      if yU>=0, th(j)(i,2)=th(j)(i,2)+w(j)*yU/nn(j); end
    end
  end

  // on-line estimation - parameters
  for j=1:nc
    thL(j)(:,t)=th(j)(:,1);
    thU(j)(:,t)=th(j)(:,2);
  end
  nu=nu+w;                              // pointer
  al=fnorm(nu);                         //   statistics and parameters
end

// RESULTS
set(scf(2),'position',[600 100 800 600])
for j=1:nc
  plot(y(1,:),y(2,:),'c.')    // data
end  
unifplot(thI,'.:r')           // initial parameters
unifplot(th)                  // estimated parameters
title 'Uniform clusters (red: initial support, blue: estimated support)'
