line source
% 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