Mercurial > hg > ltpda
diff m-toolbox/test/test_ao_polynomfit.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_polynomfit.m Wed Nov 23 19:22:13 2011 +0100 @@ -0,0 +1,199 @@ +% TEST_AO_POLYNOMFIT tests the polynomfit method of the AO class: +% - functionality +% - against lscov +% +% M Hueller and D Nicolodi 07-01-10 +% +% $Id: test_ao_polynomfit.m,v 1.2 2010/04/08 04:30:02 mauro Exp $ +% + +function test_ao_polynomfit() + %% test polynomfit vs lscov + disp(' ************** '); + disp('Example with combination of noise terms'); + disp(' '); + + fs = 1; + nsecs = 10; + n = [0 1 -2]; + + X = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm')); + N = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm')); + C = [ao(1, plist('yunits', 'm')) ao(4, plist('yunits', 'm/m')) ao(2, plist('yunits', 'm/m^(-2)'))]; + Y = C(1) * X.^0 + C(2) * X.^1 + C(3) * X.^(-2) + N; + Y.simplifyYunits; + + % Fit with UTN double/polynom_fit + fprintf('UTN POLYNOM_FIT:\n'); + try + out = polynom_fit(X.y, Y.y, n, 'nopl'); + catch err + disp(err.message); + disp('UTN method polynom_fit.m for double not available'); + end + + % Fit with Matlab lscov + fprintf('LSCOV:\n'); + x = X.y; + y = Y.y; + m = zeros(length(x), length(n)); + for k = 1:length(n) + m(:,k) = x .^ n(k); + end + + [p, stdx, mse] = lscov(m, y); + p + stdx + mse + + % Fit with LTPDA ao/polynomfit + fprintf('AO/POLYNOMFIT:\n'); + out = polynomfit(X, Y, plist('orders', n)) + + + %% More tests + %% Make fake AO from polyval + nsecs = 5; + fs = 10; + n = [0 1 -2]; + + unit_list = unit.supportedUnits; + u1 = unit(cell2mat(utils.math.randelement(unit_list,1))); + + pl1 = plist('nsecs', nsecs, 'fs', fs, ... + 'tsfcn', sprintf('t.^%d + t.^%d + t.^%d + randn(size(t))', n), ... + 'xunits', 's', 'yunits', u1); + a1 = ao(pl1); + + C = [ao(1, plist('yunits', a1.yunits)) ao(4, plist('yunits', '')) ao(2, plist('yunits', simplify(a1.yunits/(a1.yunits.^-2))))]; + N = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', u1)); + a2 = C(1) * a1.^n(1) + C(2) * a1.^n(2) + C(3) * a1.^n(3) + N; + a2.simplifyYunits; + + %% Fit with a polynomial function + + p11 = polynomfit(a1, plist('orders', n ... + )); + p12 = polynomfit(a1, plist('orders', n, ... + 'dy', 0.1)); + p22 = polynomfit(a1, plist('orders', n, ... + 'dx', 0.1, 'dy', 0.1, 'P0', zeros(size(n)))); + p13 = polynomfit(a1, plist('orders', n, ... + 'dy', 0.1*ones(size(a1.y)))); + p23 = polynomfit(a1, plist('orders', n, ... + 'dx', 0.1*ones(size(a1.x)), 'dy', 0.1*ones(size(a1.y)), 'P0', zeros(size(n)))); + p14 = polynomfit(a1, plist('orders', n, ... + 'dy', ao(0.1, plist('yunits', a1.yunits)))); + p24 = polynomfit(a1, plist('orders', n, ... + 'dx', ao(0.1, plist('yunits', a1.xunits)), 'dy', ao(0.1, plist('yunits', a1.yunits)), 'P0', zeros(size(n)))); + p15 = polynomfit(a1, plist('orders', n, ... + 'dy', ao(0.1*ones(size(a1.y)), plist('yunits', a1.yunits)))); + p25 = polynomfit(a1, plist('orders', n, ... + 'dx', ao(0.1*ones(size(a1.x)), plist('yunits', a1.xunits)), 'dy', ao(0.1*ones(size(a1.y)), plist('yunits', a1.yunits)), 'P0', zeros(size(n)))); + p26 = polynomfit(a1, plist('orders', n, ... + 'dx', 0.1, 'dy', 0.1, 'P0', ao(zeros(size(n))))); + p27 = polynomfit(a1, plist('orders', n, ... + 'dx', 0.1, 'dy', 0.1, 'P0', p11)); + + %% Compute fit: evaluating pest + + b11 = p11.eval(plist('type', 'tsdata', 'XData', a1)); + b11.removeVal(plist('value', Inf)); + b12 = p12.eval(plist('type', 'tsdata', 'XData', a1)); + b12.removeVal(plist('value', Inf)); + b22 = p22.eval(a1, plist('type', 'tsdata')); + b22.removeVal(plist('value', Inf)); + b13 = p13.eval(plist('type', 'tsdata', 'XData', a1)); + b13.removeVal(plist('value', Inf)); + b23 = p23.eval(a1, plist('type', 'tsdata')); + b23.removeVal(plist('value', Inf)); + b14 = p14.eval(plist('type', 'tsdata', 'XData', a1)); + b14.removeVal(plist('value', Inf)); + b24 = p24.eval(a1, plist('type', 'tsdata')); + b24.removeVal(plist('value', Inf)); + b15 = p15.eval(a1, plist('type', 'tsdata')); + b15.removeVal(plist('value', Inf)); + b25 = p25.eval(plist('type', 'tsdata', 'XData', a1)); + b25.removeVal(plist('value', Inf)); + b26 = p26.eval(plist('type', 'tsdata', 'XData', a1.x)); + b26.removeVal(plist('value', Inf)); + b27 = p27.eval(plist('type', 'tsdata', 'XData', a1.x)); + b27.removeVal(plist('value', Inf)); + + %% Plot fit + iplot(a1, b11, b12, b13, b14, b15, b26, b27, ... + plist('LineStyles', {'', '--'})); + + %% Remove linear trend + c = a1.removeVal(plist('value', Inf))-b11; + iplot(c) + + plot(c.hist) + + %% Reproduce from history + disp('Try rebuilding') + a_out = rebuild(c); + iplot(a_out) + plot(a_out.hist) + + %% Fit with a polynomial function + + p11 = polynomfit(a1, a2, plist('orders', n ... + )); + p12 = polynomfit(a1, a2, plist('orders', n, ... + 'dy', 0.1)); + p22 = polynomfit(a1, a2, plist('orders', n, ... + 'dx', 0.1, 'dy', 0.1, 'P0', zeros(size(n)))); + p13 = polynomfit(a1, a2, plist('orders', n, ... + 'dy', 0.1*ones(size(a1.x)))); + p23 = polynomfit(a1, a2, plist('orders', n, ... + 'dx', 0.1*ones(size(a1.x)), 'dy', 0.1*ones(size(a1.x)), 'P0', zeros(size(n)))); + p14 = polynomfit(a1, a2, plist('orders', n, ... + 'dy', ao(0.1, plist('yunits', a2.yunits)))); + p24 = polynomfit(a1, a2, plist('orders', n, ... + 'dx', ao(0.1, plist('yunits', a1.yunits)), 'dy', ao(0.1, plist('yunits', a2.yunits)), 'P0', zeros(size(n)))); + p15 = polynomfit(a1, a2, plist('orders', n, ... + 'dy', ao(0.1*ones(size(a2.y)), plist('yunits', a2.yunits)))); + p25 = polynomfit(a1, a2, plist('orders', n, ... + 'dx', ao(0.1*ones(size(a1.y)), plist('yunits', a1.yunits)), 'dy', ao(0.1*ones(size(a2.y)), plist('yunits', a2.yunits)), 'P0', zeros(size(n)))); + p26 = polynomfit(a1, a2, plist('orders', n, ... + 'dx', 0.1, 'dy', 0.1, 'P0', ao(zeros(size(n))))); + p27 = polynomfit(a1, a2, plist('orders', n, ... + 'dx', 0.1, 'dy', 0.1, 'P0', p11)); + + %% Compute fit: evaluating pest + + b11 = p11.eval(plist('type', 'xydata', 'XData', a1.y)); + b12 = p12.eval(plist('type', 'xydata', 'XData', a1)); + b22 = p22.eval(a1, plist('type', 'xydata')); + b13 = p13.eval(plist('type', 'xydata', 'XData', a1.y)); + b23 = p23.eval(a1, plist('type', 'xydata')); + b14 = p14.eval(plist('type', 'xydata', 'XData', a1.y)); + b24 = p24.eval(a1, plist('type', 'xydata')); + b15 = p15.eval(plist('type', 'xydata', 'XData', a1.y)); + b25 = p25.eval(plist('type', 'xydata', 'XData', a1.y)); + b26 = p26.eval(plist('type', 'xydata', 'XData', a1.y)); + b27 = p27.eval(plist('type', 'xydata', 'XData', a1.y)); + + % Build reference object + a12 = ao(plist('xvals', a1.y, 'yvals', a2.y, ... + 'xunits', a1.yunits, 'yunits', a2.yunits)); + %% Plot fit + iplot(a12, b11, b12, b23, b14, b25, b26, b27, ... + plist('LineStyles', {'', '--'})); + + %% Remove linear trend + c = a12-b27; + iplot(c) + + plot(c.hist) + + %% Reproduce from history + disp('Try rebuilding') + a_out = rebuild(c); + iplot(a_out) + plot(a_out.hist) + + % end + % END +end