% Example2_28.m
% modified from Example2_27.m - here use different P0 values
clear all; format short g
n=3; p=2; mu=zeros(p,1); %mu is the prior parameter estimate vector
sigma=1; sig=sigma*ones(n,1); X=[4 3; 1 2; 10 1]; y=[11 4 25]'; b=mu;
P0mag = [ 10 10^5 10^10 10^15 ];
for i = 1:length(P0mag)
    P0=P0mag(i)*eye(p,p); 
    fprintf('for P0 = %5g\n', P0mag(i));
    fprintf('  P(1,1)   P(1,2)    P(2,2)     b(1)   b(2)\n');
    P=P0; B=mu';	%P0 is V sub(beta)
    for ii=1:n;			%SEQUENTIAL
        A=P*X(ii,:)';	%eq. (2.9.8a)	%SEQUENTIAL
        Delta=sig(ii)^2+X(ii,:)*A;	%eq. (2.9.8)	%SEQUENTIAL
        K=A/Delta;	%eq. (2.9.8c)	%SEQUENTIAL
        e=y(ii)-X(ii,:)*b;	%eq. (2.9.8d)	%SEQUENTIAL
        b=b+K*e;	%eq. (2.9.8e)	%SEQUENTIAL
        P=P-K*A';	%eq. (2.9.8f)	%SEQUENTIAL
        B=[B; b'];
        if ii==1
            PP=[P(1,1) P(1,2) P(2,2)];
        else
            PP=[PP; P(1,1) P(1,2) P(2,2)];
        end
        fprintf(' %8.5g', [PP(ii,:) B(ii+1,:)]');
        fprintf('\n')
    end
end