line source
+ − % LSCOV is a wrapper for MATLAB's lscov function.
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ − %
+ − % DESCRIPTION: LSCOV is a wrapper for MATLAB's lscov function. It solves a
+ − % set of linear equations by performing a linear least-squares fit. It
+ − % solves the problem
+ − %
+ − % Y = HX
+ − %
+ − % where X are the parameters, Y the measurements, and H the linear
+ − % equations relating the two.
+ − %
+ − % CALL: X = lscov([C1 C2 ... CN], Y, pl)
+ − % X = lscov(C1,C2,C3,...,CN, Y, pl)
+ − %
+ − % INPUTS: C1...CN - AOs which represent the columns of H.
+ − % Y - AO which represents the measurement set
+ − %
+ − % Note: the length of the vectors in Ci and Y must be the same.
+ − % Note: the last input AO is taken as Y.
+ − %
+ − % pl - parameter list (see below)
+ − %
+ − % OUTPUTs: X - A pest object with fields:
+ − % y - the N fitting coefficients to y_i
+ − % dy - the parameters' standard deviations (lscov 'STDX' vector)
+ − % cov - the parameters' covariance matrix (lscov 'COV' vector)
+ − %
+ − % The procinfo field of the output PEST object is filled with the following key/value
+ − % pairs:
+ − % 'MSE' - the mean-squared errors
+ − %
+ − % <a href="matlab:utils.helper.displayMethodInfo('ao', 'lscov')">Parameters Description</a>
+ − %
+ − % VERSION: $Id: lscov.m,v 1.42 2011/04/08 08:56:17 hewitson Exp $
+ − %
+ − % EXAMPLES:
+ − %
+ − % % 1) Determine the coefficients of a linear combination of noises:
+ − %
+ − % % Make some data
+ − % fs = 10;
+ − % nsecs = 10;
+ − % B1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T'));
+ − % B2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T'));
+ − % B3 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T'));
+ − % n = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
+ − % c = [ao(1,plist('yunits','m/T')) ao(2,plist('yunits','m/T')) ao(3,plist('yunits','m T^-1'))];
+ − % y = c(1)*B1 + c(2)*B2 + c(3)*B3 + n;
+ − % y.simplifyYunits;
+ − % % Get a fit for c
+ − % p_s = lscov(B1, B2, B3, y);
+ − % % do linear combination: using lincom
+ − % yfit1 = lincom(B1, B2, B3, p_s);
+ − % yfit1.simplifyYunits;
+ − % % do linear combination: using eval
+ − % yfit2 = p_s.eval(B1, B2, B3);
+ − %
+ − % % Plot (compare data with fit)
+ − % iplot(y, yfit1, yfit2, plist('Linestyles', {'-','--'}))
+ − %
+ − % % 2) Determine the coefficients of a linear combination of noises:
+ − %
+ − % % Make some data
+ − % fs = 10;
+ − % nsecs = 10;
+ − % x1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T'));
+ − % x2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
+ − % x3 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'C'));
+ − % n = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
+ − % c = [ao(1,plist('yunits','m/T')) ao(2,plist('yunits','m/m')) ao(3,plist('yunits','m C^-1'))];
+ − % y = c(1)*x1 + c(2)*x2 + c(3)*x3 + n;
+ − % y.simplifyYunits;
+ − % % Get a fit for c
+ − % p_m = lscov(x1, x2, x3, y);
+ − % % do linear combination: using lincom
+ − % yfit1 = lincom(x1, x2, x3, p_m);
+ − % % do linear combination: using eval
+ − % pl_split = plist('times', [1 5]);
+ − % yfit2 = p_m.eval(plist('Xdata', {split(x1, pl_split), split(x2, pl_split), split(x3, pl_split)}));
+ − % % Plot (compare data with fit)
+ − % iplot(y, yfit1, yfit2, plist('Linestyles', {'-','--'}))
+ − %
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ −
+ − function varargout = lscov(varargin)
+ −
+ − % Check if this is a call for parameters
+ − if utils.helper.isinfocall(varargin{:})
+ − varargout{1} = getInfo(varargin{3});
+ − return
+ − end
+ −
+ − import utils.const.*
+ − utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
+ −
+ − % Collect input variable names
+ − in_names = cell(size(varargin));
+ − for ii = 1:nargin,in_names{ii} = inputname(ii);end
+ −
+ − % Collect all AOs and plists
+ − [A, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
+ − pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
+ −
+ − if numel(A) < 2
+ − error('### lscov needs at least 2 inputs AOs');
+ − end
+ −
+ − if nargout == 0
+ − error('### lscov can not be used as a modifier method. Please give at least one output');
+ − end
+ −
+ − % combine plists
+ − pl = parse(pl, getDefaultPlist());
+ −
+ − % collect inputs names
+ − argsname = A(1).name;
+ − for jj = 2:numel(A)
+ − argsname = [argsname ',' A(jj).name];
+ − end
+ −
+ − % Extract parameters
+ − W = find(pl, 'weights');
+ − V = find(pl, 'cov');
+ −
+ − % Build matrices for lscov
+ − C = A(1:end-1);
+ − Y = A(end);
+ −
+ − H = C(:).y;
+ − y = Y.y;
+ − if isa(W,'ao'), W = W.y; end;
+ − if isa(V,'ao'), V = V.y; end;
+ −
+ − if isempty(V)
+ − [p,stdx,mse,s] = lscov(H, y, W);
+ − else
+ − [p,stdx,mse,s] = lscov(H, y, W, V);
+ − end
+ −
+ − % prepare model, units, names
+ − model = [];
+ − for jj = 1:numel(p)
+ − names{jj} = ['C' num2str(jj)];
+ − units{jj} = Y.yunits / C(jj).yunits;
+ − xunits{jj} = C(jj).yunits;
+ − xvar{jj} = ['X' num2str(jj)];
+ − if jj == 1
+ − model = ['C' num2str(jj) '*X' num2str(jj) ' '];
+ − else
+ − model = [model ' + C' num2str(jj) '*X' num2str(jj)];
+ − end
+ − end
+ −
+ − model = smodel(plist('expression', model, ...
+ − 'params', names, ...
+ − 'values', p, ...
+ − 'xvar', xvar, ...
+ − 'xunits', xunits, ...
+ − 'yunits', Y.yunits ...
+ − ));
+ −
+ − % Build the output pest object
+ − X = pest;
+ − X.setY(p);
+ − X.setDy(stdx);
+ − X.setCov(s);
+ − X.setChi2(mse);
+ − X.setNames(names{:});
+ − X.setYunits(units{:});
+ − X.setModels(model);
+ − X.name = sprintf('lscov(%s)', argsname);
+ − X.addHistory(getInfo('None'), pl, ao_invars, [A(:).hist]);
+ − % Set procinfo object
+ − X.procinfo = plist('MSE', mse);
+ − % Propagate 'plotinfo'
+ − plotinfo = [A(:).plotinfo];
+ − if ~isempty(plotinfo)
+ − X.plotinfo = combine(plotinfo);
+ − end
+ −
+ − % Set outputs
+ − varargout{1} = X;
+ −
+ − end
+ −
+ − %--------------------------------------------------------------------------
+ − % Get Info Object
+ − %--------------------------------------------------------------------------
+ − function ii = getInfo(varargin)
+ − if nargin == 1 && strcmpi(varargin{1}, 'None')
+ − sets = {};
+ − pl = [];
+ − else
+ − sets = {'Default'};
+ − pl = getDefaultPlist;
+ − end
+ − % Build info object
+ − 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);
+ − ii.setModifier(false);
+ − ii.setArgsmin(2);
+ − end
+ −
+ − %--------------------------------------------------------------------------
+ − % Get Default Plist
+ − %--------------------------------------------------------------------------
+ − function plout = getDefaultPlist()
+ − persistent pl;
+ − if exist('pl', 'var')==0 || isempty(pl)
+ − pl = buildplist();
+ − end
+ − plout = pl;
+ − end
+ −
+ − function pl = buildplist()
+ − pl = plist();
+ −
+ − % Weights
+ − p = param({'weights','An ao containing weights for the fit.'}, paramValue.EMPTY_DOUBLE);
+ − pl.append(p);
+ −
+ − % Cov
+ − p = param({'cov','An ao containing a covariance matrix for the fit.'}, paramValue.EMPTY_DOUBLE);
+ − pl.append(p);
+ −
+ − end
+ − % END