Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/lscov.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 % LSCOV is a wrapper for MATLAB's lscov function.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: LSCOV is a wrapper for MATLAB's lscov function. It solves a
5 % set of linear equations by performing a linear least-squares fit. It
6 % solves the problem
7 %
8 % Y = HX
9 %
10 % where X are the parameters, Y the measurements, and H the linear
11 % equations relating the two.
12 %
13 % CALL: X = lscov([C1 C2 ... CN], Y, pl)
14 % X = lscov(C1,C2,C3,...,CN, Y, pl)
15 %
16 % INPUTS: C1...CN - AOs which represent the columns of H.
17 % Y - AO which represents the measurement set
18 %
19 % Note: the length of the vectors in Ci and Y must be the same.
20 % Note: the last input AO is taken as Y.
21 %
22 % pl - parameter list (see below)
23 %
24 % OUTPUTs: X - A pest object with fields:
25 % y - the N fitting coefficients to y_i
26 % dy - the parameters' standard deviations (lscov 'STDX' vector)
27 % cov - the parameters' covariance matrix (lscov 'COV' vector)
28 %
29 % The procinfo field of the output PEST object is filled with the following key/value
30 % pairs:
31 % 'MSE' - the mean-squared errors
32 %
33 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'lscov')">Parameters Description</a>
34 %
35 % VERSION: $Id: lscov.m,v 1.42 2011/04/08 08:56:17 hewitson Exp $
36 %
37 % EXAMPLES:
38 %
39 % % 1) Determine the coefficients of a linear combination of noises:
40 %
41 % % Make some data
42 % fs = 10;
43 % nsecs = 10;
44 % B1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T'));
45 % B2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T'));
46 % B3 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T'));
47 % n = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
48 % c = [ao(1,plist('yunits','m/T')) ao(2,plist('yunits','m/T')) ao(3,plist('yunits','m T^-1'))];
49 % y = c(1)*B1 + c(2)*B2 + c(3)*B3 + n;
50 % y.simplifyYunits;
51 % % Get a fit for c
52 % p_s = lscov(B1, B2, B3, y);
53 % % do linear combination: using lincom
54 % yfit1 = lincom(B1, B2, B3, p_s);
55 % yfit1.simplifyYunits;
56 % % do linear combination: using eval
57 % yfit2 = p_s.eval(B1, B2, B3);
58 %
59 % % Plot (compare data with fit)
60 % iplot(y, yfit1, yfit2, plist('Linestyles', {'-','--'}))
61 %
62 % % 2) Determine the coefficients of a linear combination of noises:
63 %
64 % % Make some data
65 % fs = 10;
66 % nsecs = 10;
67 % x1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T'));
68 % x2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
69 % x3 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'C'));
70 % n = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
71 % c = [ao(1,plist('yunits','m/T')) ao(2,plist('yunits','m/m')) ao(3,plist('yunits','m C^-1'))];
72 % y = c(1)*x1 + c(2)*x2 + c(3)*x3 + n;
73 % y.simplifyYunits;
74 % % Get a fit for c
75 % p_m = lscov(x1, x2, x3, y);
76 % % do linear combination: using lincom
77 % yfit1 = lincom(x1, x2, x3, p_m);
78 % % do linear combination: using eval
79 % pl_split = plist('times', [1 5]);
80 % yfit2 = p_m.eval(plist('Xdata', {split(x1, pl_split), split(x2, pl_split), split(x3, pl_split)}));
81 % % Plot (compare data with fit)
82 % iplot(y, yfit1, yfit2, plist('Linestyles', {'-','--'}))
83 %
84 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
85
86 function varargout = lscov(varargin)
87
88 % Check if this is a call for parameters
89 if utils.helper.isinfocall(varargin{:})
90 varargout{1} = getInfo(varargin{3});
91 return
92 end
93
94 import utils.const.*
95 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
96
97 % Collect input variable names
98 in_names = cell(size(varargin));
99 for ii = 1:nargin,in_names{ii} = inputname(ii);end
100
101 % Collect all AOs and plists
102 [A, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
103 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
104
105 if numel(A) < 2
106 error('### lscov needs at least 2 inputs AOs');
107 end
108
109 if nargout == 0
110 error('### lscov can not be used as a modifier method. Please give at least one output');
111 end
112
113 % combine plists
114 pl = parse(pl, getDefaultPlist());
115
116 % collect inputs names
117 argsname = A(1).name;
118 for jj = 2:numel(A)
119 argsname = [argsname ',' A(jj).name];
120 end
121
122 % Extract parameters
123 W = find(pl, 'weights');
124 V = find(pl, 'cov');
125
126 % Build matrices for lscov
127 C = A(1:end-1);
128 Y = A(end);
129
130 H = C(:).y;
131 y = Y.y;
132 if isa(W,'ao'), W = W.y; end;
133 if isa(V,'ao'), V = V.y; end;
134
135 if isempty(V)
136 [p,stdx,mse,s] = lscov(H, y, W);
137 else
138 [p,stdx,mse,s] = lscov(H, y, W, V);
139 end
140
141 % prepare model, units, names
142 model = [];
143 for jj = 1:numel(p)
144 names{jj} = ['C' num2str(jj)];
145 units{jj} = Y.yunits / C(jj).yunits;
146 xunits{jj} = C(jj).yunits;
147 xvar{jj} = ['X' num2str(jj)];
148 if jj == 1
149 model = ['C' num2str(jj) '*X' num2str(jj) ' '];
150 else
151 model = [model ' + C' num2str(jj) '*X' num2str(jj)];
152 end
153 end
154
155 model = smodel(plist('expression', model, ...
156 'params', names, ...
157 'values', p, ...
158 'xvar', xvar, ...
159 'xunits', xunits, ...
160 'yunits', Y.yunits ...
161 ));
162
163 % Build the output pest object
164 X = pest;
165 X.setY(p);
166 X.setDy(stdx);
167 X.setCov(s);
168 X.setChi2(mse);
169 X.setNames(names{:});
170 X.setYunits(units{:});
171 X.setModels(model);
172 X.name = sprintf('lscov(%s)', argsname);
173 X.addHistory(getInfo('None'), pl, ao_invars, [A(:).hist]);
174 % Set procinfo object
175 X.procinfo = plist('MSE', mse);
176 % Propagate 'plotinfo'
177 plotinfo = [A(:).plotinfo];
178 if ~isempty(plotinfo)
179 X.plotinfo = combine(plotinfo);
180 end
181
182 % Set outputs
183 varargout{1} = X;
184
185 end
186
187 %--------------------------------------------------------------------------
188 % Get Info Object
189 %--------------------------------------------------------------------------
190 function ii = getInfo(varargin)
191 if nargin == 1 && strcmpi(varargin{1}, 'None')
192 sets = {};
193 pl = [];
194 else
195 sets = {'Default'};
196 pl = getDefaultPlist;
197 end
198 % Build info object
199 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.op, '$Id: lscov.m,v 1.42 2011/04/08 08:56:17 hewitson Exp $', sets, pl);
200 ii.setModifier(false);
201 ii.setArgsmin(2);
202 end
203
204 %--------------------------------------------------------------------------
205 % Get Default Plist
206 %--------------------------------------------------------------------------
207 function plout = getDefaultPlist()
208 persistent pl;
209 if exist('pl', 'var')==0 || isempty(pl)
210 pl = buildplist();
211 end
212 plout = pl;
213 end
214
215 function pl = buildplist()
216 pl = plist();
217
218 % Weights
219 p = param({'weights','An ao containing weights for the fit.'}, paramValue.EMPTY_DOUBLE);
220 pl.append(p);
221
222 % Cov
223 p = param({'cov','An ao containing a covariance matrix for the fit.'}, paramValue.EMPTY_DOUBLE);
224 pl.append(p);
225
226 end
227 % END