Mercurial > hg > ltpda

comparison m-toolbox/classes/@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 % POLYNOMFIT is a polynomial fitting tool
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: POLYNOMFIT is a polynomial fitting tool based on MATLAB's
5 % lscov function. It solves an equation in the form
6 %
7 % Y = P(1) * X^N(1) + P(2) * X^N(2) + ...
8 %
9 % for the fit parameters P. It handles arbitrary powers of the input vector
10 % and uncertainties on the dependent vector Y and input vectors X.
11 % The output is a pest object where the fields are containing:
12 % Quantity % Field
13 % Fit coefficients y
14 % Uncertainties on the fit parameters
15 % (given as standard deviations) dy
16 % The reduced CHI2 of the fit chi2
17 % The covariance matrix cov
18 % The degrees of freedom of the fit dof
19 %
20 % CALL: P = polynomfit(X, Y, PL)
21 % P = polynomfit(A, PL)
22 %
23 % INPUTS: Y - dependent variable
24 % X - input variables
25 % A - data ao whose x and y fields are used in the fit
26 % PL - parameter list
27 %
28 % OUTPUT: P - a pest object with M = numel(N) fitting coefficients
29 %
30 %
31 % PARAMETERS:
32 % 'orders' - polynom orders. Eg [0,1,-2] fits to P0 + P1*x + P2./x.^2
33 % 'dy' - uncertainty on the dependent variable
34 % 'dx' - uncertainties on the input variable
35 % 'p0' - initial guess on the fit parameters used ONLY to propagate
36 % uncertainities in the input variable X to the dependent variable Y
37 %
38 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'polynomfit')">Parameters Description</a>
39 %
40 % VERSION: $Id: polynomfit.m,v 1.20 2011/04/08 08:56:11 hewitson Exp $
41 %
42 % EXAMPLES:
43 %
44 % % 1) Fit with one object input
45 %
46 % nsecs = 5;
47 % fs = 10;
48 % n = [0 1 -2];
49 % u1 = unit('mV');
50 %
51 % pl1 = plist('nsecs', nsecs, 'fs', fs, ...
52 % 'tsfcn', sprintf('t.^%d + t.^%d + t.^%d + randn(size(t))', n), ...
53 % 'xunits', 's', 'yunits', u1);
54 % a1 = ao(pl1);
55 % out1 = polynomfit(a1, plist('orders', n, 'dx', 0.1, 'dy', 0.1, 'P0', zeros(size(n))));
56 %
57 % % 2) Fit with two objects input
58 %
59 % fs = 1;
60 % nsecs = 10;
61 % n = [0 1 -2];
62 %
63 % X = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm', 'name', 'base'));
64 % N = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm', 'name', 'noise'));
65 % C = [ao(1, plist('yunits', 'm', 'name', 'C1')) ...
66 % ao(4, plist('yunits', 'm/m', 'name', 'C2')) ...
67 % ao(2, plist('yunits', 'm/m^(-2)', 'name', 'C3'))];
68 % Y = C(1) * X.^0 + C(2) * X.^1 + C(3) * X.^(-2) + N;
69 % Y.simplifyYunits;
70 % out2 = polynomfit(X, Y, plist('orders', n))
71 %
72 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
73
74 function varargout = polynomfit(varargin)
75
76 % check if this is a call for parameters
77 if utils.helper.isinfocall(varargin{:})
78 varargout{1} = getInfo(varargin{3});
79 return
80 end
81
82 % tell the system we are runing
83 import utils.const.*
84 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
85
86 % collect input variable names
87 in_names = cell(size(varargin));
88 for ii = 1:nargin,in_names{ii} = inputname(ii);end
89
90 % collect all AOs and plists
91 [aos, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
92 pli = utils.helper.collect_objects(varargin(:), 'plist', in_names);
93
94 if nargout == 0
95 error('### polynomfit can not be used as a modifier method. Please give at least one output');
96 end
97
98 % combine plists, making sure the user input is not empty
99 pli = combine(pli, plist());
100 pl = parse(pli, getDefaultPlist());
101
102 % extract arguments
103 if (length(aos) == 1)
104 % we are using x and y fields of the single ao we have
105 x = aos(1).x;
106 dx = aos(1).dx;
107 y = aos(1).y;
108 dy = aos(1).dy;
109 xunits = aos(1).xunits;
110 yunits = aos(1).yunits;
111 argsname = aos(1).name;
112 elseif (length(aos) == 2)
113 % we are using y fields of the two aos we have
114 x = aos(1).y;
115 dx = aos(1).dy;
116 y = aos(2).y;
117 dy = aos(2).dy;
118 xunits = aos(1).yunits;
119 yunits = aos(2).yunits;
120 argsname = [aos(1).name ',' aos(2).name];
121 else
122 error('### polynomfit needs one or two input AOs');
123 end
124
125 % extract plist parameters. For dx and dy we check the user input plist before
126 dy = find(pli, 'dy', dy);
127 dx = find(pli, 'dx', dx);
128 n = find(pl, 'orders');
129 p0 = find(pl, 'p0');
130
131 % vectors length
132 N = length(y);
133
134 % number of parameters
135 num_params = length(n);
136
137 % uncertainty on Y
138 if isempty(dy)
139 dy = 1;
140 end
141 if isa(dy, 'ao')
142 % check units
143 if yunits ~= dy.data.yunits
144 error('### Y and DY units are not compatible - %s %s', char(yunits), char(dy.data.yunits));
145 end
146 % extract values from AO
147 dy = dy.y;
148 end
149 if isscalar(dy)
150 % given a single value construct a vector
151 dy = ones(N, 1) * dy;
152 end
153
154 % weights
155 sigma2 = dy.^2;
156
157 % extract values for initial guess
158 if (isa(p0, 'ao') || isa(p0, 'pest'))
159 p0 = p0.y;
160 end
161
162 % uncertainty on X
163 if ~isempty(dx)
164
165 if length(p0) ~= num_params
166 error('### initial parameters guess p0 is mandatory for proper handling of X uncertainties');
167 end
168
169 if isa(dx, 'ao')
170 % check units
171 if xunits ~= dx.data.yunits
172 error('### X and DX units are not compatible - %s %s', char(xunits), char(dx.data.yunits));
173 end
174 % extract values from AO
175 dx = dx.y;
176 end
177 if isscalar(dx)
178 % given a single value construct a vector
179 dx = ones(N, 1) * dx;
180 end
181
182 % propagate X uncertainty on Y
183 dy_dx_mod = zeros(N, 1);
184 for k = 1:num_params
185 dy_dx_mod = dy_dx_mod + n(k) * p0(k) * x.^(n(k)-1);
186 end
187 sigma2x = dy_dx_mod.^2 .* dx.^2;
188
189 % add contribution to weights
190 sigma2 = sigma2 + sigma2x;
191
192 end
193
194 % construct matrix with desired powers of X
195 m = zeros(length(x), num_params);
196 for k = 1:num_params
197 m(:,k) = x .^ n(k);
198 end
199
200 % check for the presence of 1/0 terms
201 M = [];
202 X = [];
203 Y = [];
204 S = [];
205 kk = 0;
206
207 idx = isfinite(m);
208 for jj = 1:size(m,1)
209 if all(idx(jj,:))
210 kk = kk + 1;
211 M(kk,:) = m(jj,:);
212 X(kk) = x(jj);
213 Y(kk) = y(jj);
214 S(kk,:) = sigma2(jj,:);
215 end
216 end
217 m = M; clear M;
218 x = X'; clear X;
219 y = Y'; clear Y;
220 sigma2 = S; clear S;
221 N = kk;
222
223 % solve
224 [p, stdp, mse, s] = lscov(m, y, 1./sigma2);
225
226 % scale errors
227 stdp = stdp ./ sqrt(mse);
228 s = s ./ (mse);
229
230 % compute chi2
231 dof = N - length(p);
232 chi2 = sum((y - polynomeval(x, n, p)).^2 ./ sigma2) / dof;
233
234 % prepare model, units, names
235 model = [];
236 for kk = 1:length(p)
237 if kk == 1
238 model = [model 'P' num2str(kk) '*X.^(' num2str(n(kk)) ')'];
239 else
240 model = [model ' + P' num2str(kk) '*X.^(' num2str(n(kk)) ')'];
241 end
242 units(kk) = simplify(yunits/xunits.^(n(kk)));
243 names{kk} = ['P' num2str(kk)];
244 end
245 model = smodel(plist('expression', model, ...
246 'params', names, ...
247 'values', p, ...
248 'xvar', 'X', ...
249 'xunits', xunits, ...
250 'yunits', yunits ...
251 ));
252
253
254 % Build the output pest object
255 out = pest;
256 out.setY(p);
257 out.setDy(stdp);
258 out.setCov(s);
259 out.setChi2(chi2);
260 out.setDof(dof);
261 out.setNames(names{:});
262 out.setYunits(units);
263 out.setModels(model);
264 out.name = sprintf('polynomfit(%s)', argsname);
265 out.addHistory(getInfo('None'), pl, ao_invars, [aos(:).hist]);
266 % Set procinfo object
267 out.procinfo = plist('MSE', mse);
268
269 % set outputs
270 varargout{1} = out;
271
272 end
273
274 % computes polynomial combination
275 function out = polynomeval(x, n, p)
276 assert(length(p) == length(n));
277 out = zeros(size(x, 1), 1);
278 for k = 1:length(n)
279 out = out + p(k) * x.^n(k);
280 end
281 end
282
283 % get info object
284 function ii = getInfo(varargin)
285 if nargin == 1 && strcmpi(varargin{1}, 'None')
286 sets = {};
287 pl = [];
288 else
289 sets = {'Default'};
290 pl = getDefaultPlist();
291 end
292 % build info object
293 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.op, '$Id: polynomfit.m,v 1.20 2011/04/08 08:56:11 hewitson Exp $', sets, pl);
294 ii.setModifier(false);
295 ii.setArgsmin(1);
296 end
297
298 % get default plist
299 function plout = getDefaultPlist()
300 persistent pl;
301 if ~exist('pl', 'var') || isempty(pl)
302 pl = buildplist();
303 end
304 plout = pl;
305 end
306
307 function pl = buildplist()
308 pl = plist();
309
310 % orders
311 p = param({'orders', 'Polynom orders.'}, [0]);
312 pl.append(p);
313
314 % default plist for linear fitting
315 pl.append(plist.LINEAR_FIT_PLIST);
316
317 end