Mercurial > hg > ltpda
diff m-toolbox/test/test_ao_bilinfit.m @ 0:f0afece42f48
Import.
| author | Daniele Nicolodi <nicolodi@science.unitn.it> |
|---|---|
| date | Wed, 23 Nov 2011 19:22:13 +0100 |
| parents | |
| children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/m-toolbox/test/test_ao_bilinfit.m Wed Nov 23 19:22:13 2011 +0100 @@ -0,0 +1,208 @@ +% Test ao/bilinfit for: +% - functionality +% - against lscov +% +% M Hueller and D Nicolodi 07-01-10 +% +% $Id: test_ao_bilinfit.m,v 1.5 2010/05/07 14:04:42 nicolodi Exp $ +% +function test_ao_bilinfit() + + + %% test bilinfit vs lscov + disp(' ************** '); + disp('Example with combination of noise terms'); + disp(' '); + fs = 10; + nsecs = 10; + x1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T')); + x2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T')); + n = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm')); + c = [ao(1,plist('yunits','m/T')) ao(2,plist('yunits','m/T'))]; + y = c(1)*x1 + c(2)*x2 + n; + y.simplifyYunits; + + % Get a fit for c, using lscov (no constant term) + disp('Using ao/lscov, no constant term'); + disp(' '); + p_lscov = lscov(x1, x2, y); + disp(' '); + p_lscov.y + disp(' '); + + % do linear combination + yfit_lscov = lincom(x1, x2, p_lscov); + + % Get a fit for c, using bilinfit + disp('Using ao/bilinfit'); + disp(' '); + p_b = bilinfit(x1, x2, y); + disp(' '); + p_b.y + disp(' '); + disp(' ************** '); + + % do linear combination: use direct sum + yfit_b1 = simplifyYunits(x1 * find(p_b, 'P1') + x2 * find(p_b, 'P2') + find(p_b, 'P3')); + % do linear combination: use eval + yfit_b2 = p_b.eval(plist('Xdata', {x1, x2})); + + % Plot (compare data with fit) + iplot(y, yfit_lscov, yfit_b1, yfit_b2, plist('Linestyles', {'all','-'})) + + %% test with uncertainties + fprintf('\n\n\n'); + + fs = 1; + nsecs = 50; + x1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'C')); + x2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm')); + n = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'N')); + c = ao(plist('xvals', x1.x, 'yvals', ones(size(x1.y)), 'type', 'tsdata')); + b = [ao(1,plist('yunits','N/C')) ao(2,plist('yunits','N/m'))]; + y = b(1)*x1 + b(2)*x2 + n; + y.simplifyYunits(); + + p_m = lscov(x1, x2, c, y, plist('weights', ao(plist('vals', 1./(ones(size(x1.x))))))); + fprintf('AO LSCOV: Fit results: A = %f +/- %f, B = %f +/- %f, Y0 = %f +/- %f\n', ... + p_m.y(1), p_m.dy(1)./sqrt(find(p_m.procinfo, 'mse')), ... + p_m.y(2), p_m.dy(2)./sqrt(find(p_m.procinfo, 'mse')), ... + p_m.y(3), p_m.dy(3)./sqrt(find(p_m.procinfo, 'mse'))); + + % do linear combination: using linear combination + yfit1 = simplifyYunits(x1 * find(p_m, 'C1') + x2 * find(p_m, 'C2') + find(p_m, 'C3')); + + % do linear combination: using eval + yfit2 = p_m.eval(plist('Xdata',{x1, x2, c})); + + % Plot (compare data with fit) + iplot(y, yfit1, yfit2, plist('Linestyles', {'-','--'})) + + fprintf('UTN BILINFIT:\n'); + try + [p, perr] = bilinfit(y.y, x1.y, x2.y, ones(size(x1.x))); + catch err + disp(err.message); + disp('UTN method bilinfit.m for double not available'); + end + [x,stdx,mse] = lscov([c.y, x1.y, x2.y], y.y, 1./(ones(size(x1.x)))); + fprintf('M LSCOV: Fit results: A = %f +/- %f, B = %f +/- %f, Y0 = %f +/- %f\n', ... + x(2), stdx(2)*sqrt(1/mse), ... + x(3), stdx(3)*sqrt(1/mse), ... + x(1), stdx(1)*sqrt(1/mse)); + + p_b = bilinfit(x1, x2, y, plist('dy', ao(plist('vals', ones(size(x1.x)), 'yunits', 'N')))); + fprintf('BILINFIT: Fit results: A = %f +/- %f, B = %f +/- %f, Y0 = %f +/- %f\n', ... + p_b.y(1), p_b.dy(1), ... + p_b.y(2), p_b.dy(2), ... + p_b.y(3), p_b.dy(3)); + + %% very simple test + x1 = ao(plist('xvals', [1 2 3], 'yvals', [1 2 3])); + x2 = ao(plist('xvals', [1 2 3], 'yvals', [6 4 7])); + y = x1 + x2; + c = ao(plist('xvals', x1.x, 'yvals', ones(size(x1.y)), 'type', 'tsdata')); + + % Fit with UTN double/bilinfit + fprintf('UTN BILINFIT:\n'); + try + [p, perr] = bilinfit(y.y, x1.y, x2.y, ones(size(x1.x))); + catch err + disp(err.message); + disp('UTN method bilinfit.m for double not available'); + end + + % Fit with Matlab lscov + [x,stdx,mse] = lscov([c.y, x1.y, x2.y], y.y, 1./(ones(size(x1.x)))); + fprintf('M LSCOV: Fit results: A = %f +/- %f, B = %f +/- %f, Y0 = %f +/- %f\n', ... + x(2), stdx(2)*sqrt(1/mse), ... + x(3), stdx(3)*sqrt(1/mse), ... + x(1), stdx(1)*sqrt(1/mse)); + + % Fit with LTPDA ao/bilinfit + [x] = bilinfit(x1, x2, y, plist('dy', ao(plist('vals', ones(size(x1.x)))))); + fprintf('BILINFIT: Fit results: A = %f +/- %f, B = %f +/- %f, Y0 = %f +/- %f\n', ... + x.y(1), x.dy(1), ... + x.y(2), x.dy(2), ... + x.y(3), x.dy(3)); + + %% More tests about format of inputs + + %% Make fake AO + nsecs = 100; + fs = 10; + + unit_list = unit.supportedUnits; + u1 = unit(cell2mat(utils.math.randelement(unit_list,1))); + u2 = unit(cell2mat(utils.math.randelement(unit_list,1))); + u3 = unit(cell2mat(utils.math.randelement(unit_list,1))); + + pl1 = plist('nsecs', nsecs, 'fs', fs, ... + 'tsfcn', 'randn(size(t))', ... + 'yunits', u1); + + pl2 = plist('nsecs', nsecs, 'fs', fs, ... + 'tsfcn', 'randn(size(t))', ... + 'yunits', u2); + + pl3 = plist('nsecs', nsecs, 'fs', fs, ... + 'tsfcn', 'randn(size(t))', ... + 'yunits', u3); + + x1 = ao(pl1); + x2 = ao(pl2); + x3 = ao(pl3); + + n = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm')); + c = [ao(1,plist('yunits', unit('m')/u1)) ao(2,plist('yunits', unit('m')/u2)) ao(3,plist('yunits', unit('m')/u3))]; + y = c(1)*x1 + c(2)*x2 + c(3)*x3 + n; + y.simplifyYunits; + + %% Fit with bilinfit + + p11 = bilinfit(x1, x2, x3, y, plist(... + )); + p12 = bilinfit(x1, x2, x3, y, plist(... + 'dy', 0.1 ... + )); + p22 = bilinfit(x1, x2, x3, y, plist(... + 'dx', [0.1 0.1 0.1], ... + 'dy', 0.1, ... + 'P0', [0 0 0 0])); + p13 = bilinfit(x1, x2, x3, y, plist(... + 'dy', ao(0.1, plist('yunits', y.yunits)) ... + )); + p23 = bilinfit(x1, x2, x3, y, plist(... + 'dx', [ao(0.1, plist('yunits', x1.yunits)) ... + ao(0.1, plist('yunits', x2.yunits)) ... + ao(0.1, plist('yunits', x3.yunits))], ... + 'dy', ao(0.1, plist('yunits', y.yunits)), ... + 'P0', [0 0 0 0])); + p14 = bilinfit(x1, x2, x3, y, plist(... + 'dy', ao(0.1*ones(size(x1.y)), plist('yunits', y.yunits)) ... + )); + p24 = bilinfit(x1, x2, x3, y, plist(... + 'dx', [ao(0.1*ones(size(x1.x)), plist('yunits', x1.yunits)) ... + ao(0.1*ones(size(x2.x)), plist('yunits', x2.yunits)) ... + ao(0.1*ones(size(x3.x)), plist('yunits', x3.yunits))], ... + 'dy', ao(0.1*ones(size(x1.y)), plist('yunits', y.yunits)), ... + 'P0', [0 0 0 0])); + p25 = bilinfit(x1, x2, x3, y, plist(... + 'dx', [0.1 0.1 0.1], ... + 'dy', 0.1, ... + 'P0', ao([0 0 0 0]))); + p26 = bilinfit(x1, x2, x3, y, plist(... + 'dx', [0.1 0.1 0.1], ... + 'dy', 0.1, ... + 'P0', p11)); + + %% Compute fit: evaluating pest + + b11 = p11.eval(plist('XData', {x1, x2, x3})); + b12 = p12.eval(x1, x2, x3); + b22 = p22.eval(plist('XData', {x1.y, x2.y, x3.y})); + b23 = p23.eval(plist('XData', [x1 x2 x3])); + b24 = p24.eval(x1, x2, x3); + b25 = p25.eval(plist('XData', {x1, x2, x3})); + b26 = p26.eval([x1 x2 x3]); +end