--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/straight_line_fit/straight_line_fit.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,148 @@
+% 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