% Tutorial: Bayesian Linear Regression

close all;
% set seed for random number generator
rand('seed',1);

% generate data.
numExamples = 10;
xGrid = (-1:0.001:1)';
x = randsample(xGrid,numExamples);
% polynom: y = x^3 + 0.7x^2-0.5x + 0.2 plus
y = x.^3+0.7*x.^2-0.5*x+0.2;
% add some noise
y = y + 0.05*randn(length(x),1);

% plot data
scatter(x,y,'filled');
axis([-1 1 0 1]);
hold on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TASK 1: Plot original polynom
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
XGrid = [ones(size(xGrid)) xGrid xGrid.^2 xGrid.^3];
w = [0.2 -0.5 0.7 1];       % added
origY = XGrid * w';         % added
orig = plot(xGrid,origY);   % added
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% generate polynomial features to train linear regression
X = [ones(size(x)) x x.^2 x.^3];

% posterior of bayesian linear regression (see lecture): normal
% distribution with parameter thetaStar and A^-1
I = eye(size(x,2));
sigma = 0.05;
sigma_p = 1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TASK 2: Compute MAP-model and covariance-matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
thetaStar = (X'*X+sigma^2/sigma_p^2*I)\X'*y;    % slide 43
A = sigma^-2*(X'*X)+sigma_p^-2*I;               % slide 43
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% draw sample from posterior
numSamples = 3;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TASK 3: Sample polynoms from posterior
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
W = mvnrnd(thetaStar,inv(A),numSamples);        % slide 53 ff
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TASK 4: Plot sampled polynoms
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
XGrid = [ones(size(xGrid)) xGrid xGrid.^2 xGrid.^3];    % added
sampled = plot(xGrid,XGrid*W');                         % added
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% TASK 5: Use different parameter for numExamples

