Mercurial > hg > ltpda

comparison m-toolbox/test/test_ao_polynomfit.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % TEST_AO_POLYNOMFIT tests the polynomfit method of the AO class:
2 % - functionality
3 % - against lscov
4 %
5 % M Hueller and D Nicolodi 07-01-10
6 %
7 % $Id: test_ao_polynomfit.m,v 1.2 2010/04/08 04:30:02 mauro Exp $
8 %
9
10 function test_ao_polynomfit()
11 %% test polynomfit vs lscov
12 disp(' ************** ');
13 disp('Example with combination of noise terms');
14 disp(' ');
15
16 fs = 1;
17 nsecs = 10;
18 n = [0 1 -2];
19
20 X = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
21 N = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
22 C = [ao(1, plist('yunits', 'm')) ao(4, plist('yunits', 'm/m')) ao(2, plist('yunits', 'm/m^(-2)'))];
23 Y = C(1) * X.^0 + C(2) * X.^1 + C(3) * X.^(-2) + N;
24 Y.simplifyYunits;
25
26 % Fit with UTN double/polynom_fit
27 fprintf('UTN POLYNOM_FIT:\n');
28 try
29 out = polynom_fit(X.y, Y.y, n, 'nopl');
30 catch err
31 disp(err.message);
32 disp('UTN method polynom_fit.m for double not available');
33 end
34
35 % Fit with Matlab lscov
36 fprintf('LSCOV:\n');
37 x = X.y;
38 y = Y.y;
39 m = zeros(length(x), length(n));
40 for k = 1:length(n)
41 m(:,k) = x .^ n(k);
42 end
43
44 [p, stdx, mse] = lscov(m, y);
45 p
46 stdx
47 mse
48
49 % Fit with LTPDA ao/polynomfit
50 fprintf('AO/POLYNOMFIT:\n');
51 out = polynomfit(X, Y, plist('orders', n))
52
53
54 %% More tests
55 %% Make fake AO from polyval
56 nsecs = 5;
57 fs = 10;
58 n = [0 1 -2];
59
60 unit_list = unit.supportedUnits;
61 u1 = unit(cell2mat(utils.math.randelement(unit_list,1)));
62
63 pl1 = plist('nsecs', nsecs, 'fs', fs, ...
64 'tsfcn', sprintf('t.^%d + t.^%d + t.^%d + randn(size(t))', n), ...
65 'xunits', 's', 'yunits', u1);
66 a1 = ao(pl1);
67
68 C = [ao(1, plist('yunits', a1.yunits)) ao(4, plist('yunits', '')) ao(2, plist('yunits', simplify(a1.yunits/(a1.yunits.^-2))))];
69 N = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', u1));
70 a2 = C(1) * a1.^n(1) + C(2) * a1.^n(2) + C(3) * a1.^n(3) + N;
71 a2.simplifyYunits;
72
73 %% Fit with a polynomial function
74
75 p11 = polynomfit(a1, plist('orders', n ...
76 ));
77 p12 = polynomfit(a1, plist('orders', n, ...
78 'dy', 0.1));
79 p22 = polynomfit(a1, plist('orders', n, ...
80 'dx', 0.1, 'dy', 0.1, 'P0', zeros(size(n))));
81 p13 = polynomfit(a1, plist('orders', n, ...
82 'dy', 0.1*ones(size(a1.y))));
83 p23 = polynomfit(a1, plist('orders', n, ...
84 'dx', 0.1*ones(size(a1.x)), 'dy', 0.1*ones(size(a1.y)), 'P0', zeros(size(n))));
85 p14 = polynomfit(a1, plist('orders', n, ...
86 'dy', ao(0.1, plist('yunits', a1.yunits))));
87 p24 = polynomfit(a1, plist('orders', n, ...
88 'dx', ao(0.1, plist('yunits', a1.xunits)), 'dy', ao(0.1, plist('yunits', a1.yunits)), 'P0', zeros(size(n))));
89 p15 = polynomfit(a1, plist('orders', n, ...
90 'dy', ao(0.1*ones(size(a1.y)), plist('yunits', a1.yunits))));
91 p25 = polynomfit(a1, plist('orders', n, ...
92 '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))));
93 p26 = polynomfit(a1, plist('orders', n, ...
94 'dx', 0.1, 'dy', 0.1, 'P0', ao(zeros(size(n)))));
95 p27 = polynomfit(a1, plist('orders', n, ...
96 'dx', 0.1, 'dy', 0.1, 'P0', p11));
97
98 %% Compute fit: evaluating pest
99
100 b11 = p11.eval(plist('type', 'tsdata', 'XData', a1));
101 b11.removeVal(plist('value', Inf));
102 b12 = p12.eval(plist('type', 'tsdata', 'XData', a1));
103 b12.removeVal(plist('value', Inf));
104 b22 = p22.eval(a1, plist('type', 'tsdata'));
105 b22.removeVal(plist('value', Inf));
106 b13 = p13.eval(plist('type', 'tsdata', 'XData', a1));
107 b13.removeVal(plist('value', Inf));
108 b23 = p23.eval(a1, plist('type', 'tsdata'));
109 b23.removeVal(plist('value', Inf));
110 b14 = p14.eval(plist('type', 'tsdata', 'XData', a1));
111 b14.removeVal(plist('value', Inf));
112 b24 = p24.eval(a1, plist('type', 'tsdata'));
113 b24.removeVal(plist('value', Inf));
114 b15 = p15.eval(a1, plist('type', 'tsdata'));
115 b15.removeVal(plist('value', Inf));
116 b25 = p25.eval(plist('type', 'tsdata', 'XData', a1));
117 b25.removeVal(plist('value', Inf));
118 b26 = p26.eval(plist('type', 'tsdata', 'XData', a1.x));
119 b26.removeVal(plist('value', Inf));
120 b27 = p27.eval(plist('type', 'tsdata', 'XData', a1.x));
121 b27.removeVal(plist('value', Inf));
122
123 %% Plot fit
124 iplot(a1, b11, b12, b13, b14, b15, b26, b27, ...
125 plist('LineStyles', {'', '--'}));
126
127 %% Remove linear trend
128 c = a1.removeVal(plist('value', Inf))-b11;
129 iplot(c)
130
131 plot(c.hist)
132
133 %% Reproduce from history
134 disp('Try rebuilding')
135 a_out = rebuild(c);
136 iplot(a_out)
137 plot(a_out.hist)
138
139 %% Fit with a polynomial function
140
141 p11 = polynomfit(a1, a2, plist('orders', n ...
142 ));
143 p12 = polynomfit(a1, a2, plist('orders', n, ...
144 'dy', 0.1));
145 p22 = polynomfit(a1, a2, plist('orders', n, ...
146 'dx', 0.1, 'dy', 0.1, 'P0', zeros(size(n))));
147 p13 = polynomfit(a1, a2, plist('orders', n, ...
148 'dy', 0.1*ones(size(a1.x))));
149 p23 = polynomfit(a1, a2, plist('orders', n, ...
150 'dx', 0.1*ones(size(a1.x)), 'dy', 0.1*ones(size(a1.x)), 'P0', zeros(size(n))));
151 p14 = polynomfit(a1, a2, plist('orders', n, ...
152 'dy', ao(0.1, plist('yunits', a2.yunits))));
153 p24 = polynomfit(a1, a2, plist('orders', n, ...
154 'dx', ao(0.1, plist('yunits', a1.yunits)), 'dy', ao(0.1, plist('yunits', a2.yunits)), 'P0', zeros(size(n))));
155 p15 = polynomfit(a1, a2, plist('orders', n, ...
156 'dy', ao(0.1*ones(size(a2.y)), plist('yunits', a2.yunits))));
157 p25 = polynomfit(a1, a2, plist('orders', n, ...
158 '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))));
159 p26 = polynomfit(a1, a2, plist('orders', n, ...
160 'dx', 0.1, 'dy', 0.1, 'P0', ao(zeros(size(n)))));
161 p27 = polynomfit(a1, a2, plist('orders', n, ...
162 'dx', 0.1, 'dy', 0.1, 'P0', p11));
163
164 %% Compute fit: evaluating pest
165
166 b11 = p11.eval(plist('type', 'xydata', 'XData', a1.y));
167 b12 = p12.eval(plist('type', 'xydata', 'XData', a1));
168 b22 = p22.eval(a1, plist('type', 'xydata'));
169 b13 = p13.eval(plist('type', 'xydata', 'XData', a1.y));
170 b23 = p23.eval(a1, plist('type', 'xydata'));
171 b14 = p14.eval(plist('type', 'xydata', 'XData', a1.y));
172 b24 = p24.eval(a1, plist('type', 'xydata'));
173 b15 = p15.eval(plist('type', 'xydata', 'XData', a1.y));
174 b25 = p25.eval(plist('type', 'xydata', 'XData', a1.y));
175 b26 = p26.eval(plist('type', 'xydata', 'XData', a1.y));
176 b27 = p27.eval(plist('type', 'xydata', 'XData', a1.y));
177
178 % Build reference object
179 a12 = ao(plist('xvals', a1.y, 'yvals', a2.y, ...
180 'xunits', a1.yunits, 'yunits', a2.yunits));
181 %% Plot fit
182 iplot(a12, b11, b12, b23, b14, b25, b26, b27, ...
183 plist('LineStyles', {'', '--'}));
184
185 %% Remove linear trend
186 c = a12-b27;
187 iplot(c)
188
189 plot(c.hist)
190
191 %% Reproduce from history
192 disp('Try rebuilding')
193 a_out = rebuild(c);
194 iplot(a_out)
195 plot(a_out.hist)
196
197 % end
198 % END
199 end