Mercurial > hg > ltpda

comparison m-toolbox/test/test_ao_bilinfit.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/bilinfit for:
2 % - functionality
3 % - against lscov
4 %
5 % M Hueller and D Nicolodi 07-01-10
6 %
7 % $Id: test_ao_bilinfit.m,v 1.5 2010/05/07 14:04:42 nicolodi Exp $
8 %
9 function test_ao_bilinfit()
10
11
12 %% test bilinfit vs lscov
13 disp(' ************** ');
14 disp('Example with combination of noise terms');
15 disp(' ');
16 fs = 10;
17 nsecs = 10;
18 x1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T'));
19 x2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T'));
20 n = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
21 c = [ao(1,plist('yunits','m/T')) ao(2,plist('yunits','m/T'))];
22 y = c(1)*x1 + c(2)*x2 + n;
23 y.simplifyYunits;
24
25 % Get a fit for c, using lscov (no constant term)
26 disp('Using ao/lscov, no constant term');
27 disp(' ');
28 p_lscov = lscov(x1, x2, y);
29 disp(' ');
30 p_lscov.y
31 disp(' ');
32
33 % do linear combination
34 yfit_lscov = lincom(x1, x2, p_lscov);
35
36 % Get a fit for c, using bilinfit
37 disp('Using ao/bilinfit');
38 disp(' ');
39 p_b = bilinfit(x1, x2, y);
40 disp(' ');
41 p_b.y
42 disp(' ');
43 disp(' ************** ');
44
45 % do linear combination: use direct sum
46 yfit_b1 = simplifyYunits(x1 * find(p_b, 'P1') + x2 * find(p_b, 'P2') + find(p_b, 'P3'));
47 % do linear combination: use eval
48 yfit_b2 = p_b.eval(plist('Xdata', {x1, x2}));
49
50 % Plot (compare data with fit)
51 iplot(y, yfit_lscov, yfit_b1, yfit_b2, plist('Linestyles', {'all','-'}))
52
53 %% test with uncertainties
54 fprintf('\n\n\n');
55
56 fs = 1;
57 nsecs = 50;
58 x1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'C'));
59 x2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
60 n = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'N'));
61 c = ao(plist('xvals', x1.x, 'yvals', ones(size(x1.y)), 'type', 'tsdata'));
62 b = [ao(1,plist('yunits','N/C')) ao(2,plist('yunits','N/m'))];
63 y = b(1)*x1 + b(2)*x2 + n;
64 y.simplifyYunits();
65
66 p_m = lscov(x1, x2, c, y, plist('weights', ao(plist('vals', 1./(ones(size(x1.x)))))));
67 fprintf('AO LSCOV: Fit results: A = %f +/- %f, B = %f +/- %f, Y0 = %f +/- %f\n', ...
68 p_m.y(1), p_m.dy(1)./sqrt(find(p_m.procinfo, 'mse')), ...
69 p_m.y(2), p_m.dy(2)./sqrt(find(p_m.procinfo, 'mse')), ...
70 p_m.y(3), p_m.dy(3)./sqrt(find(p_m.procinfo, 'mse')));
71
72 % do linear combination: using linear combination
73 yfit1 = simplifyYunits(x1 * find(p_m, 'C1') + x2 * find(p_m, 'C2') + find(p_m, 'C3'));
74
75 % do linear combination: using eval
76 yfit2 = p_m.eval(plist('Xdata',{x1, x2, c}));
77
78 % Plot (compare data with fit)
79 iplot(y, yfit1, yfit2, plist('Linestyles', {'-','--'}))
80
81 fprintf('UTN BILINFIT:\n');
82 try
83 [p, perr] = bilinfit(y.y, x1.y, x2.y, ones(size(x1.x)));
84 catch err
85 disp(err.message);
86 disp('UTN method bilinfit.m for double not available');
87 end
88 [x,stdx,mse] = lscov([c.y, x1.y, x2.y], y.y, 1./(ones(size(x1.x))));
89 fprintf('M LSCOV: Fit results: A = %f +/- %f, B = %f +/- %f, Y0 = %f +/- %f\n', ...
90 x(2), stdx(2)*sqrt(1/mse), ...
91 x(3), stdx(3)*sqrt(1/mse), ...
92 x(1), stdx(1)*sqrt(1/mse));
93
94 p_b = bilinfit(x1, x2, y, plist('dy', ao(plist('vals', ones(size(x1.x)), 'yunits', 'N'))));
95 fprintf('BILINFIT: Fit results: A = %f +/- %f, B = %f +/- %f, Y0 = %f +/- %f\n', ...
96 p_b.y(1), p_b.dy(1), ...
97 p_b.y(2), p_b.dy(2), ...
98 p_b.y(3), p_b.dy(3));
99
100 %% very simple test
101 x1 = ao(plist('xvals', [1 2 3], 'yvals', [1 2 3]));
102 x2 = ao(plist('xvals', [1 2 3], 'yvals', [6 4 7]));
103 y = x1 + x2;
104 c = ao(plist('xvals', x1.x, 'yvals', ones(size(x1.y)), 'type', 'tsdata'));
105
106 % Fit with UTN double/bilinfit
107 fprintf('UTN BILINFIT:\n');
108 try
109 [p, perr] = bilinfit(y.y, x1.y, x2.y, ones(size(x1.x)));
110 catch err
111 disp(err.message);
112 disp('UTN method bilinfit.m for double not available');
113 end
114
115 % Fit with Matlab lscov
116 [x,stdx,mse] = lscov([c.y, x1.y, x2.y], y.y, 1./(ones(size(x1.x))));
117 fprintf('M LSCOV: Fit results: A = %f +/- %f, B = %f +/- %f, Y0 = %f +/- %f\n', ...
118 x(2), stdx(2)*sqrt(1/mse), ...
119 x(3), stdx(3)*sqrt(1/mse), ...
120 x(1), stdx(1)*sqrt(1/mse));
121
122 % Fit with LTPDA ao/bilinfit
123 [x] = bilinfit(x1, x2, y, plist('dy', ao(plist('vals', ones(size(x1.x))))));
124 fprintf('BILINFIT: Fit results: A = %f +/- %f, B = %f +/- %f, Y0 = %f +/- %f\n', ...
125 x.y(1), x.dy(1), ...
126 x.y(2), x.dy(2), ...
127 x.y(3), x.dy(3));
128
129 %% More tests about format of inputs
130
131 %% Make fake AO
132 nsecs = 100;
133 fs = 10;
134
135 unit_list = unit.supportedUnits;
136 u1 = unit(cell2mat(utils.math.randelement(unit_list,1)));
137 u2 = unit(cell2mat(utils.math.randelement(unit_list,1)));
138 u3 = unit(cell2mat(utils.math.randelement(unit_list,1)));
139
140 pl1 = plist('nsecs', nsecs, 'fs', fs, ...
141 'tsfcn', 'randn(size(t))', ...
142 'yunits', u1);
143
144 pl2 = plist('nsecs', nsecs, 'fs', fs, ...
145 'tsfcn', 'randn(size(t))', ...
146 'yunits', u2);
147
148 pl3 = plist('nsecs', nsecs, 'fs', fs, ...
149 'tsfcn', 'randn(size(t))', ...
150 'yunits', u3);
151
152 x1 = ao(pl1);
153 x2 = ao(pl2);
154 x3 = ao(pl3);
155
156 n = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
157 c = [ao(1,plist('yunits', unit('m')/u1)) ao(2,plist('yunits', unit('m')/u2)) ao(3,plist('yunits', unit('m')/u3))];
158 y = c(1)*x1 + c(2)*x2 + c(3)*x3 + n;
159 y.simplifyYunits;
160
161 %% Fit with bilinfit
162
163 p11 = bilinfit(x1, x2, x3, y, plist(...
164 ));
165 p12 = bilinfit(x1, x2, x3, y, plist(...
166 'dy', 0.1 ...
167 ));
168 p22 = bilinfit(x1, x2, x3, y, plist(...
169 'dx', [0.1 0.1 0.1], ...
170 'dy', 0.1, ...
171 'P0', [0 0 0 0]));
172 p13 = bilinfit(x1, x2, x3, y, plist(...
173 'dy', ao(0.1, plist('yunits', y.yunits)) ...
174 ));
175 p23 = bilinfit(x1, x2, x3, y, plist(...
176 'dx', [ao(0.1, plist('yunits', x1.yunits)) ...
177 ao(0.1, plist('yunits', x2.yunits)) ...
178 ao(0.1, plist('yunits', x3.yunits))], ...
179 'dy', ao(0.1, plist('yunits', y.yunits)), ...
180 'P0', [0 0 0 0]));
181 p14 = bilinfit(x1, x2, x3, y, plist(...
182 'dy', ao(0.1*ones(size(x1.y)), plist('yunits', y.yunits)) ...
183 ));
184 p24 = bilinfit(x1, x2, x3, y, plist(...
185 'dx', [ao(0.1*ones(size(x1.x)), plist('yunits', x1.yunits)) ...
186 ao(0.1*ones(size(x2.x)), plist('yunits', x2.yunits)) ...
187 ao(0.1*ones(size(x3.x)), plist('yunits', x3.yunits))], ...
188 'dy', ao(0.1*ones(size(x1.y)), plist('yunits', y.yunits)), ...
189 'P0', [0 0 0 0]));
190 p25 = bilinfit(x1, x2, x3, y, plist(...
191 'dx', [0.1 0.1 0.1], ...
192 'dy', 0.1, ...
193 'P0', ao([0 0 0 0])));
194 p26 = bilinfit(x1, x2, x3, y, plist(...
195 'dx', [0.1 0.1 0.1], ...
196 'dy', 0.1, ...
197 'P0', p11));
198
199 %% Compute fit: evaluating pest
200
201 b11 = p11.eval(plist('XData', {x1, x2, x3}));
202 b12 = p12.eval(x1, x2, x3);
203 b22 = p22.eval(plist('XData', {x1.y, x2.y, x3.y}));
204 b23 = p23.eval(plist('XData', [x1 x2 x3]));
205 b24 = p24.eval(x1, x2, x3);
206 b25 = p25.eval(plist('XData', {x1, x2, x3}));
207 b26 = p26.eval([x1 x2 x3]);
208 end