% PURPOSE : To estimate the weights of an MLP with inputs x
%           and outputs y generated by a time-varying, noisy  
%           and chaotic map. The MLP is trained with both the EKF and 
%           EKF with evidence maximisation.                           
             
% AUTHOR  : Nando de Freitas - Thanks for the acknowledgement :-)
% DATE    : 08-09-98

clear;
echo off;

% INITIALISATION AND PARAMETERS:
% =============================
l=300;              % Number of time samples;
N=l-1;
t=1:1:l-1;
s1=10;              % Number of neurons in the first layer.
s2=1;               % Number of neurons in the second layer.
initVar= 1;         % Variance of prior for weights.
KalmanR = 1e-4;     % Kalman filter measurement noise covariance hyperparameter;
KalmanQ = 0;        % Kalman filter process noise covariance hyperparameter;
KalmanP = 100;      % Kalman filter initial weights covariance hyperparameter;
window  = 3;        % Length of moving window to estimate the time covariance. 


% GENERATE INPUT-OUTPUT DATA:
% ==========================
b1=3.5*ones(1,l/2);
b2=3.7*ones(1,l/4);
b3=3.1*ones(1,l/4);
b=[b1 b2 b3];
sn=0.3*ones(l,1);
for i=1:l-1,
 sn(i+1,1) = b(1,i)*sn(i,1)*(1-sn(i,1));
end;
x=sn(1:l-1,1)';                   % Input.
y=sn(2:l,1)';                     % Output.
y = y + KalmanR*randn(size(y));   % Add noise to output.
trials=1;                         % Number of simulations.


for i=1:trials,
  % PERFORM EXTENDED KALMAN FILTERING TO TRAIN MLP:
  % ==============================================
  fprintf('\n')
  fprintf('Training the MLP with EKF')
  fprintf('\n')
  flops(0);
  [p,theta,thetaR,PR,OutputVariance] = mlpekf(x,y,s1,s2,KalmanR,KalmanQ,KalmanP,initVar,N);
  durationekf=flops;
  errorekf=norm(p-y);


  % PERFORM SEQUENTIAL EKF WITH Q UPDATE:
  % ====================================
  fprintf('\n')
  fprintf('Training the MLP with EKFMAX')
  fprintf('\n')
  fprintf('\n')
  flops(0);
  [p,theta,thetaR,PR,OutputVariance,qplot] = mlpekfQ(x,y,s1,s2,KalmanR,KalmanQ,KalmanP,initVar,window,N);
  durationekf2=flops;
  errorekf2=norm(p-y);

  fprintf('EKF error    = %6.2f',errorekf)
  fprintf('\n')
  fprintf('EKFMAX error = %6.2f ',errorekf2)
  fprintf('\n')
  fprintf('EKF duration    = %8.2f flops.\n',durationekf)
  fprintf('EKFMAX duration = %8.2f flops.\n',durationekf2) 
  fprintf('\n')

  % PLOT RESULTS:
  % ============
											
  figure(1)
  clf;
  line=-.1:.01:1.2;
  line2=1e-4:1e-4:3e-4;
  subplot(311)
  plot(1:length(y),y,'k:',1:length(y),p,'r');
  ylabel('Prediction','fontsize',15)
  axis([0 300 0.1 1.2]); 
  hold on
  plot(150*ones(size(line)),line,'k');
  plot(225*ones(size(line)),line,'k');
  axis([0 300 0.1 1.2]); 
  subplot(312)
  plot(1:length(OutputVariance),OutputVariance,'m');
  ylabel('Output variance','fontsize',15);
  axis([0 300 0.0001 0.00022]);
  hold on;
  plot(150*ones(size(line2)),line2,'k');
  plot(225*ones(size(line2)),line2,'k');
  axis([0 300 0.0001 0.00022]);
  subplot(313)
  plot(1:length(qplot),qplot,'b');
  ylabel('Noise parameter','fontsize',15);
  axis([0 300 -0.001 0.012]); 
  hold on;
  plot(150*ones(size(line)),line,'k');
  plot(225*ones(size(line)),line,'k');
  axis([0 300 -0.001 0.012])
end;