% 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; endfunction 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