% Example 2_32.m Gauss iteration
b = [ 2 1]';  % initial guesses
psi = [ 1 1 1; 1 2 2; 1 2 3 ];
U = [ 0 0; 0 4];
mu = b;
W = inv(psi);
% data
t = [-1 0 1]';
Y = [1 3 8.5]';
n = length(Y);
p = 2;
% initialize based on the initial values of b
X = [ ones(n,1) 2*(1+b(2)*t).*t ];
eta = b(1) + (1+b(2)*t).^2;
S = (Y-eta)'*W*(Y-eta)+(mu-b)'*U*(mu-b);
fprintf(' k+1   b_1        b_2    deltab_1/b1    deltab_2/b2      S\n')
fprintf(' %2d  %8.6f  %8.6f                             %8.6g\n',...
           0, b', S );
for i = 1:4
    Pinv = X'*W*X + U;
    delta_b = Pinv \ (X'*W*(Y - eta)+ U*(mu - b ));
    b_errs = (delta_b ./ b );
    b = b + delta_b;
    % now update with the new b
    X = [ ones(n,1) 2*(1+b(2)*t).*t ];
    eta = b(1) + (1+b(2)*t).^2;
    S = (Y-eta)'*W*(Y-eta)+(mu-b)'*U*(mu-b);
    fprintf(' %2d  %8.6f  %8.6f  %11.6g  %11.6g  %11.7g\n',...
        i, b', b_errs', S );

end