%
% Linear regression example
%

%%
% Generate the data
%
%  yk = theta1 + theta2 tk

  theta1 = 1;
  theta2 = 0.5;
  sd = 0.1;
  
  t  = (0:0.01:1)';
  ry = theta1 + theta2 * t;
  y  = ry + sd * randn(size(ry));
  
  h = plot(t,ry,'-',t,y,'.');
  set(h,'LineWidth',2);
  set(h,'MarkerSize',15);
  
%%
% Batch linear regression
%
  m0 = [0;0];
  P0 = eye(2);
  H  = [ones(size(y)) t];
  
  mb = inv(inv(P0) + 1/sd^2*H'*H) * (1/sd^2*H'*y + inv(P0)*m0)
  Pb = inv(inv(P0) + 1/sd^2*H'*H)
  
  ey = mb(1) + mb(2) * t;
  h = plot(t,ry,'-',t,y,'.',t,ey,'--');
  set(h,'LineWidth',2);
  set(h,'MarkerSize',15);

%%
% Recursive linear regression
%
  m = [0;0];
  P = eye(2);
  
  for k=1:length(y)
      Hk = [1 t(k)];
      m = inv(inv(P) + 1/sd^2*Hk'*Hk) * (1/sd^2*Hk'*y(k) + inv(P)*m)
      P = inv(inv(P) + 1/sd^2*Hk'*Hk)
      m
      pause
  end
  
%%
% Recursive linear regression II
%
  m = [0;0];
  P = eye(2);
  
  for k=1:length(y)
      Hk = [1 t(k)];
      S = Hk * P * Hk' + sd^2;
      K = P * Hk' / S;
      m = m + K * (y(k) - Hk * m);
      P = P - K * S * K';
  end
  m