Mercurial > hg > ltpda

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