% this is the testing Matlab file of pooling for Experiment Two 
% described in the paper "Can MC Supplement
% Adaptive AR Predictors" published % in CMP96. 
% Subfiles includes esti1.m; codec.m;
% code1.m; code2.m


clear all
global m T theta y my lambda p t ty np
my=3;                               % number of the states of Y
T=1000;                             % overall simulation length
m=3;                                % memory length
lambda=0.99;                        % forgetting factor
np=3;                               % number of predictors
fbar=1/my;                          % pre-prior pdf
y=zeros(1,T);
rand('seed',39);
counter=50;
thetabad=zeros(my,counter*my^m);
errorm=0;
error0=0;
model=500;
Xbar=eye(1,m);                    % expectation of regressor
Xbar0=Xbar;
sigma=10;
P=sigma*eye(m,m);                      % P=(X-EX)*(X-EX)'
P0=P;
lambda=0.99;
nu=m;                                  % forgetting time


 error_yell=zeros(1,model);    % ideal model error
 error_red=zeros(1,model);   % pooled model error
 error_blue=zeros(1,model);
 error_green=zeros(1,model);
 error_white=zeros(1,model);

for mod=1:model
v0=(1/np)*ones(np,1);                % pooling weight
v=v0*ones(1,T);

% the transition probability matrix-theta
theta0=rand(m,my^m);
for j=1:my^m
  thetas(j)=0;
  for i=1:m
  thetas(j)=thetas(j)+theta0(i,j);
  end
end
for i=1:m
for j=1:my^m
theta(i,j)=theta0(i,j)/thetas(j);
end
end

% initial states

for i=1:m
  y(i)=floor(rand*my+1);
end

% generating data y(t)

for t=m+1:T
  ty=codec(y);
  rn=rand;
  pom=0;
  for f=1:my
     if (pom<rn)&(pom+theta(f,ty)>rn)
       y(t)=f;               
     end
     pom=pom+theta(f,ty);
  end  
end

  sita=ones(my,np*my^m);                     % estimated value of theta  
  tf1=sita;
  tf2=sita;
  epsi=1;
  nc=epsi*ones(my,np*my^m);
  n1=nc;
  n2=nc;
  prec=ones(my,1);
  pred1=ones(my,np);
  pred2=pred1;
  pictc=zeros(1,T); 
  pict1=pictc;
  pict2=pictc;

for t=m+1:T
    for i=1:m
     X(i)=y(t-i);
    end
     nu=lambda*(nu+1);
     Xbar=lambda*Xbar+(X-Xbar)/(nu)+(1-lambda)*Xbar0;     
     P=lambda*P+((X-Xbar)'*(X-Xbar)-P)/(nu)+(1-lambda)*P0;
     R=chol(P);
  % prediction 
  ty=codec(y); 
  p=1;
  [sita,nc]=esti1(nc,sita);             %  correct predictor 
  prec=sita(:,ty);  
  pictc(t)=expect(prec);
 for p=1:np
 ty=filt1(X,p,R,m);
%    ty=code1(y);
  % updating pooling weight v
v(p,t)=(tf1(y(t),ty+(p-1)*my^m)*v(p,t))^lambda*v0(p)^(1-lambda);
  
    [tf1,n1]=esti1(n1,tf1);          % approximate predictors
    pred1(:,p)=tf1(:,ty);
    ty=code2(y);
    [tf2,n2]=esti1(n2,tf2);
    pred2(:,p)=tf2(:,ty);
  end
 v(:,t)=v(:,t)/norm(v(:,t),1);
  pict11(t)=expect(pred1(:,[1]));
  pict12(t)=expect(pred1(:,[2]));
  pict13(t)=expect(pred1(:,[3]));
   
   % pooling

   sum1=0;
   sum2=0;
   for i=1:my
       f01(i)=1;
       f02(i)=1;
        for j=1:np
        f01(i)=f01(i)*pred1(i,j)^(v(j,t));
        f02(i)=f02(i)*pred2(i,j)^(v(j,t));
        end
        sum1=sum1+f01(i);
        sum2=sum2+f02(i); 
   end
   f01=f01/sum1;
   f02=f02/sum2;
   pict1(t)=expect(f01);
   pict2(t)=expect(f02);
   
end
frac=4;                        % initial part omitted
if(frac<=T/(T-2)),frac=T/(T-2);end
tt=(T/frac):1:T;
picc=pictc(tt);
pic1=pict1(tt);
pic2=pict2(tt);
pic11=pict11(tt);
pic12=pict12(tt);
pic13=pict13(tt);
yy=y(tt);
error_yell(1,mod)=sqrt((yy-picc)*(yy-picc)'*frac/(frac-1)/T);
error_red(1,mod)=sqrt((yy-pic1)*(yy-pic1)'*frac/(frac-1)/T);
error_blue(1,mod)=sqrt((yy-pic12)*(yy-pic12)'*frac/(frac-1)/T);
error_green(1,mod)=sqrt((yy-pic11)*(yy-pic11)'*frac/(frac-1)/T);
error_white(1,mod)=sqrt((yy-pic13)*(yy-pic13)'*frac/(frac-1)/T);

if error_red(1,mod)>min([error_blue(1,mod) error_green(1,mod) error_white(1,mod)])
error0=error0+1;
if error0<counter, thetabad(:,(error0-1)*my^m+1:error0*my^m)=theta;end;  
end
end
quanlity=error_red./error_yell;
quanlity1=error_blue./error_yell;
error0

qq=min([error_blue;error_green;error_white]);
quanlity2=qq./error_red;