% eda05_04
% illustrate least squares fitting of a line
% when the data have unequal variance
% this version uses F and H

% set of vector of x's
N=51;
N2=(N-1)/2;
Dx = 2;
xmin = 0.0;
xmax = xmin+Dx*(N-1);
x = xmin+Dx*[0:N-1]';

% prepare some test data hat consists of
% two sections, that differ slightly in
% interrcept and slope
a = 1.0;
b = 2.0;
d = a + b*x;
Da = 1;
Db = 0.3;
d(N2:N) = d(N2:N) + Da + Db*x(N2:N);

% prepare the variance, high for thr first section
% low for the second section
sd = zeros(N,1);
sd(1:N2-1) = 100;
sd(N2:N) = 5;

% set up the data kernel
G=zeros(N,2);
G(:,1)=ones(N,1);
G(:,2)=x;

% set up F and f
F=zeros(N,2);
for i=[1:N]
    F(i,:)=G(i,:)./sd(i);
end
f=d./sd;

% solve for estimated model parameters

mest = (F'*F)\(F'*f);
% predict data
dpre = mest(1) + mest(2) .* x;

figure(1);
clf;
set(gca,'LineWidth',2);
hold on;
axis( [0, 100, -500, 500] );
plot( x, d, 'ko', 'LineWidth', 2 );

% plot 2-sigma error bars
for p = [1:N]
    plot( [x(p), x(p)], [d(p)-2*sd(p), d(p)+2*sd(p)], 'k-', 'LineWidth', 2 );
    hold on;
    end
plot(x,dpre,'k-','LineWidth',2);

% now swap the variances
sd = zeros(N,1);
sd(1:N2-1) = 5;
sd(N2:N) = 100;
sd2i = sd.^(-2);
Cdi = diag(sd2i);

% set up the data kernel
G=zeros(N,2);
G(:,1)=ones(N,1);
G(:,2)=x;

% set up F and f
F=zeros(N,2);
for i=[1:N]
    F(i,:)=G(i,:)./sd(i);
end
f=d./sd;

% solve for estimated model parameters
mest = (F'*F)\(F'*f);
% predict data
dpre = mest(1) + mest(2) .* x;

figure(2);
clf;
set(gca,'LineWidth',2);
hold on;
axis( [0, 100, -500, 500] );
plot( x, d, 'ko', 'LineWidth', 2 );

% plot 2-sigma error bars
for p = [1:N]
    plot( [x(p), x(p)], [d(p)-2*sd(p), d(p)+2*sd(p)], 'k-', 'LineWidth', 2 );
    hold on;
    end
plot(x,dpre,'k-','LineWidth',2);
hold on;




