MATLAB Examples

Contents

function bayesian_regression(X,Y,small_sigma_squared,eta_sqaured)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Name - bayesian_regression
% Creation Date - 3rd Nov 2014
% Author: Soumya Banerjee
% Website: https://sites.google.com/site/neelsoumya/
%
% Description:
%   Function to do Bayesian regression
%   inspired by video on bayesian linear regression
%   https://www.youtube.com/watch?v=qz2U8coNwV4
%   by mathematicalmonk on youtube
%
% Input:
%       X - matrix of predictors
%       Y - vector of responses
%       small_sigma_squared - standard deviation^2 (variance) for covariance matrix for Y
%       eta_sqaured - standard deviation^2 (variance) for covariance matrix for beta (regressors)
%
% Output:
%       1) Vector of inferred regressors/parameters
%       2) Histograms of inferred regressors/parameters
%       3) Monte Carlo trace plots
%
% Assumptions -
%
% Example usage:
%   X = randn(100,5)
%   r = [0;2;0;-3;0] % only two nonzero coefficients
%   Y = X*r + randn(100,1)*.1 % small added noise
%   small_sigma_squared = 0.01
%   eta_sqaured = 0.01
%   bayesian_regression(X,Y,small_sigma_squared,eta_sqaured)
%
% License - BSD
%
% Acknowledgements -
%           Dedicated to my wife Joyeeta Ghose and my mother Kalyani
%           Banerjee
%
% Change History -
%                   3rd  Nov 2014  - Creation by Soumya Banerjee
%                   20th Nov 2014  - Modified by Soumya Banerjee
%                                       no burn-in required; hence taken
%                                       out (thanks to suggestion by
%                                       Alireza Kashani)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


tic;

do Bayesian regression

beta is vector of regressors P(beta|D) ~ N(beta| mu, lambda) mu = lambda * X{transpose} * big_sigma^(-1) * Y lambda = (X{transpose} * inv(sigma) * X + inv(big_omega))^{-1} $P(\beta|D) \sim N(\beta| \mu, \lambda)$ $\mu = \Lambda * X^{T} * \Sigma^{-1} * Y$ $\Lambda = (X^{T} * \Sigma^{-1} * X + \Omega^{-1})^{-1}$

iNumMeasurements = size(X,1);
iNumRegressors   = size(X,2);

big_sigma = small_sigma_squared * eye(iNumMeasurements); %a = 0.01 * eye(5);
big_omega = eta_sqaured * eye(iNumRegressors); %b = 0.01 * eye(5);

disp('covariance matrix and mean vector of posterior distribution')
lambda = inv(X' * inv(big_sigma) * X + inv(big_omega))
mu     = lambda * X' * inv(big_sigma) * Y

% mvnrnd(MU,SIGMA)
% mvnrnd(mu,lambda)
covariance matrix and mean vector of posterior distribution

lambda =

   1.0e-03 *

    0.0930   -0.0041   -0.0036   -0.0036   -0.0106
   -0.0041    0.1179   -0.0003    0.0049    0.0176
   -0.0036   -0.0003    0.0847    0.0016   -0.0006
   -0.0036    0.0049    0.0016    0.1027    0.0018
   -0.0106    0.0176   -0.0006    0.0018    0.1031


mu =

    0.0076
    1.9786
    0.0035
   -2.9569
   -0.0123

draw samples from posterior

iNumIter = 10000; % number of samples
for iCount=1:iNumIter
    w_vector_array(iCount,:) = [mvnrnd(mu,lambda)];
end

get mean from posterior

size(w_vector_array)
disp('inferred parameter vector (mean)')
[mean(w_vector_array(1:end,1)) mean(w_vector_array(1:end,2)) ...
    mean(w_vector_array(1:end,3))  mean(w_vector_array(1:end,4)) ...
    mean(w_vector_array(1:end,5)) ]
ans =

       10000           5

inferred parameter vector (mean)

ans =

    0.0077    1.9785    0.0035   -2.9569   -0.0122

plot histograms

if iNumRegressors == 5
    iNumBins = 100;
    figID = figure;
    subplot(2,3,1)
    hist(w_vector_array(1:end,1),iNumBins)
    hold on
    subplot(2,3,2)
    hist(w_vector_array(1:end,2),iNumBins)
    hold on
    subplot(2,3,3)
    hist(w_vector_array(1:end,3),iNumBins)
    hold on
    subplot(2,3,4)
    hist(w_vector_array(1:end,4),iNumBins)
    hold on
    subplot(2,3,[5 6])
    hist(w_vector_array(1:end,5),iNumBins)

    hold off

save plots

    print(figID, '-djpeg', sprintf('bayesregression_parameters_hist%s.jpg', date));

plot Monte Carlo trace plots

    figID_2 = figure
    plot(w_vector_array(1:end,1))
    xlabel('Monte Carlo samples'); ylabel('Posterior of one regressor parameter')
    print(figID_2, '-djpeg', sprintf('mcmctrace_%s.jpg', date));
figID_2 = 

  Figure (2) with properties:

      Number: 2
        Name: ''
       Color: [0.9400 0.9400 0.9400]
    Position: [360 278 560 420]
       Units: 'pixels'

  Use GET to show all properties

end

toc;
Elapsed time is 3.993003 seconds.