Mercurial > hg > ltpda

comparison m-toolbox/classes/@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 % BILINFIT is a linear fitting tool
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: BILINFIT linear fitting tool based on MATLAB's lscov
5 % function. It solves an equation in the form
6 %
7 % Y = X(1) * P(1) + X(2) * P(2) + ... + P(N+1)
8 %
9 % for the fit parameters P. It handles an arbitrary number of input vectors
10 % and uncertainties on the dependent vector Y and input vectors X(1..N).
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 = bilinfit(X1, X2, .., XN, Y, PL)
21 %
22 % INPUTS: Y - dependent variable
23 % X(1..N) - input variables
24 % PL - parameter list
25 %
26 % OUTPUT: P - a pest object with the N+1 elements
27 %
28 %
29 % PARAMETERS:
30 % 'dy' - uncertainty on the dependent variable
31 % 'dx' - uncertainties on the input variables
32 % 'p0' - initial guess on the fit parameters to propagate uncertainities
33 % in the input variables X(1..N) to the dependent variable Y
34 %
35 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'bilinfit')">Parameters Description</a>
36 %
37 % VERSION: $Id: bilinfit.m,v 1.20 2011/04/08 08:56:11 hewitson Exp $
38 %
39 % EXAMPLES:
40 %
41 % % 1) Determine the coefficients of a linear combination of noises:
42 %
43 % % Make some data
44 % fs = 10;
45 % nsecs = 10;
46 % x1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
47 % x2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
48 % n = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
49 % c = [ao(1,plist('yunits','m/m')) ao(2,plist('yunits','m/m'))];
50 % y = c(1)*x1 + c(2)*x2 + n;
51 % y.simplifyYunits;
52 %
53 % % Get a fit for the c coefficients and a constant term
54 % p = bilinfit(x1, x2, y)
55 %
56 % % Do linear combination: using eval
57 % pl_split = plist('times', [1 5]);
58 % yfit = eval(p, split(x1, pl_split), split(x2, pl_split));
59 %
60 % % Plot (compare data with fit)
61 % iplot(y, yfit, plist('Linestyles', {'-','--'}))
62 %
63 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
64
65 function varargout = bilinfit(varargin)
66
67 % check if this is a call for parameters
68 if utils.helper.isinfocall(varargin{:})
69 varargout{1} = getInfo(varargin{3});
70 return
71 end
72
73 % tell the system we are runing
74 import utils.const.*
75 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
76
77 % collect input variable names
78 in_names = cell(size(varargin));
79 for ii = 1:nargin,in_names{ii} = inputname(ii);end
80
81 % collect all AOs and plists
82 [aos, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
83 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
84
85 if nargout == 0
86 error('### bilinfit can not be used as a modifier method. Please give at least one output');
87 end
88
89 if numel(aos) < 2
90 error('### bilinfit needs at least 2 inputs AOs');
91 end
92
93 % combine plists
94 pl = parse(pl, getDefaultPlist());
95
96 % extract parameters
97 dy = find(pl, 'dy');
98 dx = find(pl, 'dx');
99 p0 = find(pl, 'p0');
100
101 % collect inputs
102 Y = aos(end);
103 X = aos(1:end-1);
104
105 % collect inputs names
106 argsname = aos(1).name;
107 for jj = 2:numel(aos)
108 argsname = [argsname ',' aos(jj).name];
109 end
110
111 % get data from AOs
112 x = X(:).y;
113 y = Y.y;
114
115 % vectors length
116 N = length(y);
117
118 % uncertainty on Y
119 if isempty(dy)
120 dy = 1;
121 end
122 if isa(dy, 'ao')
123 % check units
124 if Y.yunits ~= dy.yunits
125 error('### Y and DY units are not compatible - %s %s', char(Y.yunits), char(dy.yunits));
126 end
127 % extract values from AO
128 dy = dy.y;
129 end
130 if isscalar(dy)
131 % given a single value construct a vector
132 dy = ones(N, 1) * dy;
133 end
134
135 % squares
136 sigma2 = dy.^2;
137 sigma2y_rms = sqrt(sum(sigma2)/N);
138
139 % extract values for initial guess
140 if (isa(p0, 'ao') || isa(p0, 'pest'))
141 p0 = p0.y;
142 end
143
144 % uncertainty on X
145 if ~isempty(dx)
146
147 for k = 1:length(dx)
148 dxi = dx(k);
149
150 if ~isempty(dxi)
151 if isa(dxi, 'ao')
152 % check units
153 if X(k).yunits ~= dxi.yunits
154 error('### X and DX units are not compatible - %s %s', char(X.yunits), char(dxi.yunits));
155 end
156 % extract values from AO
157 dxi = dxi.y;
158 end
159 if isscalar(dxi)
160 % given a single value construct a vector
161 dxi = ones(N, 1) * dxi;
162 end
163
164 % squares
165 sigma2xi = dxi.^2;
166
167 % if A0 guess are not given
168 if isempty(p0(k))
169 % set it to obtain equal error contribution to the Y error
170 sigma2xi_rms = sqrt(sum(sigma2xi)/N);
171 p0(k) = sigma2y_rms/sigma2xi_rms;
172 end
173
174 % add contribution to weights
175 sigma2 = sigma2 + sigma2xi * p0(k)^2;
176 end
177
178 end
179 end
180
181 % constant term
182 c = ones(N, 1);
183
184 % build matrix
185 m = [x c];
186
187 % solve
188 [p, stdx, mse, s] = lscov(m, y, 1./sigma2);
189
190 % scale errors and covariance matrix
191 stdp = stdx ./ sqrt(mse);
192 s = s ./ mse;
193
194 % compute chi2
195 dof = N - length(p);
196 chi2 = sum((y - lincom(m, p)).^2 ./ sigma2) / dof;
197
198 % prepare model, units, names
199 model = [];
200 for kk = 1:length(p)
201 switch kk
202 case 1
203 units(kk) = simplify(Y.yunits/X(kk).yunits);
204 model = ['P' num2str(kk) '*X' num2str(kk)];
205 xvar{kk} = ['X' num2str(kk)];
206 xunits{kk} = X(kk).yunits;
207 case length(p)
208 units(kk) = Y.yunits;
209 model = [model ' + P' num2str(kk)];
210 otherwise
211 units(kk) = simplify(Y.yunits/X(kk).yunits);
212 model = [model ' + P' num2str(kk) '*X' num2str(kk)];
213 xvar{kk} = ['X' num2str(kk)];
214 xunits{kk} = X(kk).yunits;
215 end
216 names{kk} = ['P' num2str(kk)];
217
218 end
219 model = smodel(plist('expression', model, ...
220 'params', names, ...
221 'values', p, ...
222 'xvar', xvar, ...
223 'xunits', xunits, ...
224 'yunits', Y.yunits));
225
226
227 % build the output pest object
228 out = pest;
229 out.setY(p);
230 out.setDy(stdp);
231 out.setCov(s);
232 out.setChi2(chi2);
233 out.setDof(dof);
234 out.setNames(names{:});
235 out.setYunits(units);
236 out.setModels(model);
237 out.name = sprintf('bilinfit(%s)', argsname);
238 out.addHistory(getInfo('None'), pl, ao_invars, [aos(:).hist]);
239 % set procinfo object
240 out.procinfo = plist('MSE', mse);
241
242 % set outputs
243 varargout{1} = out;
244
245 end
246
247 % computes linear combination
248 function out = lincom(x, p)
249 assert(size(x, 2) == length(p));
250 out = zeros(size(x, 1), 1);
251 for k = 1:length(p)
252 out = out + x(:,k) * p(k);
253 end
254 end
255
256 % get info object
257 function ii = getInfo(varargin)
258 if nargin == 1 && strcmpi(varargin{1}, 'None')
259 sets = {};
260 pl = [];
261 else
262 sets = {'Default'};
263 pl = getDefaultPlist();
264 end
265 % build info object
266 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.op, '$Id: bilinfit.m,v 1.20 2011/04/08 08:56:11 hewitson Exp $', sets, pl);
267 ii.setModifier(false);
268 ii.setArgsmin(2);
269 end
270
271 % get default plist
272
273 function plout = getDefaultPlist()
274 persistent pl;
275 if ~exist('pl', 'var') || isempty(pl)
276 pl = buildplist();
277 end
278 plout = pl;
279 end
280
281 function pl = buildplist()
282
283 % default plist for linear fitting
284 pl = plist.MULTILINEAR_FIT_PLIST;
285
286 end