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