% Make some data that is a straight-line then fit with curvefit and
% lscov.
% 
% 
mc

m = 1.3;
c = 4;
x = 1:30;
ca = ao(plist('dtype', 'cdata', 'yvals', c.*ones(size(x))));
xa = ao(plist('dtype', 'cdata', 'yvals', x));
  
n = 1*randn(size(x));
sn = 2;
na = sn.*ao(plist('dtype', 'cdata', 'yvals', n));
ya = m.*xa + ca + na;

% Fit with straightLineFit

fy = straightLineFit(ya);
sigma = find(fy.procinfo, 'Sigma');
params = find(sigma.procinfo, 'P')
Pstd = find(sigma.procinfo, 'STD')
Pcov = find(sigma.procinfo, 'cov')
  

% theoretical errors
vn = var(na.y);
N = length(x);
delta = N.*sum(x.^2) - sum(x).^2;

em = sqrt(N.*vn / delta)
ec = sqrt(vn .* sum(x.^2./delta))

params
Pstd = 100.*find(sigma.procinfo, 'STD')./params


%% curvefit needs an xy data
xy = ya.convert('to xydata');
pl = plist('Function', 'P(1).*Xdata + P(2)', 'P0', [1 0]);
b = curvefit(xy, pl);


%% With weights of 1/sigma^2

m = 1.3;
c = 4;
x = 1:30;
ca = ao(plist('dtype', 'cdata', 'yvals', c.*ones(size(x))));
xa = ao(plist('dtype', 'cdata', 'yvals', x));
  
n = 1*randn(size(x));
sn = 1;
na = sn.*ao(plist('dtype', 'cdata', 'yvals', n));
ya = m.*xa + ca + na;

% Fit with straightLineFit
vn = var(na.y);
W  = 1./vn*ones(size(x));
fy = straightLineFit(ya, plist('Weights', W));
sigma = find(fy.procinfo, 'Sigma');
params = find(sigma.procinfo, 'P')
Pstd = 100*find(sigma.procinfo, 'STD')./params
Pcov = find(sigma.procinfo, 'cov')

% theoretcical errors
y = ya.y.';
delta = (sum(W).*sum(W.*(x.^2))  - sum(W.*x).^2 );
tc = (sum(W.*(x.^2)) .* sum(W.*y) - sum(W.*x).*sum(W.*x.*y))./delta

find(sigma.procinfo, 'STD')
em = sqrt(sum(W)./delta)
ec = sqrt(sum(W.*(x.^2)) ./ delta)

params
Pstd = 100.*find(sigma.procinfo, 'STD')./params


%% curvefit needs an xy data
xy = ya.convert('to xydata');
pl = plist('Function', 'P(1).*Xdata + P(2)', 'P0', [1 0], 'Yerr', vn);
b = curvefit(xy, pl);

%%

M = [];
C = [];
dM = [];
dC = [];

noise = logspace(-2,1,100);

m = 1.3;
c = 4;
x = 1:30;
ca = ao(plist('dtype', 'cdata', 'yvals', c.*ones(size(x))));
xa = ao(plist('dtype', 'cdata', 'yvals', x));

for k=1:length(noise)
  
  k
  
  n = noise(k)*randn(size(x));
  na = ao(plist('dtype', 'cdata', 'yvals', n));
  ya = m.*xa + ca + na;
  
  % Fit with straightLineFit
  
  fy = straightLineFit(ya);
  sigma = find(fy.procinfo, 'Sigma');
  params = find(sigma.procinfo, 'P');
  std = find(sigma.procinfo, 'STD');
  cov = find(sigma.procinfo, 'cov');

  M = [M params(1)];
  C = [C params(2)];
  dM = [dM std(1)];
  dC = [dC std(2)];
  
  
end
  
iplot(fy,ya, plist('YerrL', {sigma,[]}, 'YerrU', {sigma,[]}))

%% Plot
  
figure
subplot(2,1,1)
shadedplot(noise, M-dM, M+dM, [0.6 0.6 0.6], [1 1 1]);
set(gca, 'XScale', 'log')
hold on
semilogx(noise, M);
xscale('log')
xlabel('Noise level (sigma)');
ylabel('gradient');

subplot(2,1,2)
shadedplot(noise, C-dC, C+dC, [0.6 0.6 0.6], [1 1 1]);
set(gca, 'XScale', 'log')
hold on
semilogx(noise, C);
xscale('log')
xlabel('Noise level (sigma)');
ylabel('Intercept');

  
% END