%************************************************************
% File:         bwr.m
% Date:         August 21, 2006
% Author:       Jo-Anne Ting
% Description:
%  Bayesian weighted regression (BWR) for automatic
%  outlier detection on a globally linear function
%
% Inputs:
%
%  (where N = number of data samples,
%         d = dimensionality of input data)
%
%  X               : training input (Nxd matrix)
%  Y               : training output (Nx1 matrix)
%  n_iter          : max number of EM iterations
%  tol             : tolerance for termination of 
%                    algorithm (proportional change in 
%                    likelihood)
%
% Outputs:
%
% A data structure "s" containing the final converged
% values of:
%
%  s.beta_sigma20 : initial covariance for beta
%  s.beta_sigma2  : covariance of beta
%  s.beta_beta0   : initial mean of regression estimate 
%  s.beta         : regression estimate
%  s.w_a0         : parameter a for distribution of w
%  s.w_b0         : parameter b for distribution of w
%  s.w            : posterior mean of weight
%  s.sigma2       : variance of output
%
%************************************************************
function [s] = bwr(X, Y, n_iter, tol)

% Augment input data with offset dimension of 1s
X1 = [X ones(length(X),1)];
[N,d] = size(X1);

% Setup of priors
%-------------------------------------------------------
% Regression coefficient 
s.beta_sigma20 = 1.e3*eye(d); 
s.beta_beta0 = zeros(d,1);       
s.beta = s.beta_beta0;

% Start with weights having a prior weight of 1  
s.w_a0 = 1;
s.w_b0 = ones(N,1);      
s.w = s.w_a0 ./ s.w_b0;   

% Initial noise variance
s.sigma2  = 1;

% EM-Iterations
%-------------------------------------------------------
for i=1:n_iter,
  
  % Re-estimate the regression coefficients
  %-----------------------------------------------------
  sum_wxx = ( X1.*(s.w*ones(1,d)) )'*X1;
  sum_wyx = ( X1.*(s.w*ones(1,d)) )'*Y;

  s.beta_sigma2 = inv(sum_wxx./s.sigma2 ...
                   + inv(s.beta_sigma20));
  s.beta = s.beta_sigma2 * (sum_wyx./s.sigma2 ...
            + inv(s.beta_sigma20)*s.beta_beta0);
  
  % Calculate weights  
  %-----------------------------------------------------
  s.w_a = s.w_a0 + 0.5;
  s.w_b = s.w_b0 + (0.5/s.sigma2) ...
          .* ((Y - X1*s.beta).^2) ...
          + 0.5.*sum((X1*s.beta_sigma2).*X1, 2);
  s.w = s.w_a./s.w_b;

  % Update s.sigma2 
  %-----------------------------------------------------
  s.sigma2 = sum(s.w .*((Y - X1*s.beta).^2)) / N +...
             sum(s.w .*(sum((X1*s.beta_sigma2).*X1,2))) / N;
  
  % TO DO: Monitor the incomplete data likelihood to
  % determine when to stop iterating

end  % for i=1:n_iter
