Mercurial > hg > ltpda
view m-toolbox/test/straight_line_fit/straight_line_fit.m @ 17:7afc99ec5f04 database-connection-manager
Update ao_model_retrieve_in_timespan
author | Daniele Nicolodi <nicolodi@science.unitn.it> |
---|---|
date | Mon, 05 Dec 2011 16:20:06 +0100 |
parents | f0afece42f48 |
children |
line wrap: on
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