% Example2_26.m
% model equation is one parameter y = X*beta
% estimate beta sequentially starting with prior information:
%    mu = b_0 = 5 with variance V = 0.25
%
% given data:
X = [ 1 2 3 4 10 ]';
y = [ 3 12 16 17 47 ]';
s = [ 1 2 2 3 4 ]';
n = length( y );  % the number of measurements

mu = 5;
V = 0.25;

% basic algorithm - used in function at end of this file
% b = zeros( n );
% N = mu/V;  D = 1/V;
% for k = 1:n
%     N = N + X(k)*y(k)*s(k)^2;
%     D = D + X(k)^2*s(k)^2;
%     b(k) = N/D;
% end

% OLS case - set V = Inf; s = ones(n,1)
% Let V = "large"
fprintf('OLS:\n')
b_OLS = b_MAP_1p( X, y, ones(n,1), mu, 1e10 );
fprintf('b_OLS:  ')
fprintf(' %6.4f', b_OLS);
fprintf('\n')

% ML case - set V = Inf; but use the actual s.d.
% Let V = "large"
fprintf('ML:\n')
b_ML = b_MAP_1p( X, y, s, mu, 1e10 );
fprintf('b_ML:  ')
fprintf(' %6.4f', b_ML);
fprintf('\n')

% MAP case - use all the given information
fprintf('MAP:\n')
b_MAP = b_MAP_1p( X, y, s, mu, V );
fprintf('b_MAP:  ')
fprintf(' %6.4f', b_MAP);
fprintf('\n')


function [b] = b_MAP_1p( X, y, s, mu, V )
    n = length( y );
    b = zeros( n,1 );
    N = mu/V;  D = 1/V;
    for k = 1:n
        N = N + X(k)*y(k)/s(k)^2;
        D = D + X(k)^2/s(k)^2;
        b(k) = N/D;
    end
end
