Mercurial > hg > ltpda
view m-toolbox/classes/@ao/bilinfit.m @ 16:91f21a0aab35 database-connection-manager
Update utils.jquery
* * *
Update utils.jmysql.getsinfo
author | Daniele Nicolodi <nicolodi@science.unitn.it> |
---|---|
date | Mon, 05 Dec 2011 16:20:06 +0100 |
parents | f0afece42f48 |
children |
line wrap: on
line source
% BILINFIT is a linear fitting tool %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % DESCRIPTION: BILINFIT linear fitting tool based on MATLAB's lscov % function. It solves an equation in the form % % Y = X(1) * P(1) + X(2) * P(2) + ... + P(N+1) % % for the fit parameters P. It handles an arbitrary number of input vectors % and uncertainties on the dependent vector Y and input vectors X(1..N). % The output is a pest object where the fields are containing: % Quantity % Field % Fit coefficients y % Uncertainties on the fit parameters % (given as standard deviations) dy % The reduced CHI2 of the fit chi2 % The covariance matrix cov % The degrees of freedom of the fit dof % % CALL: P = bilinfit(X1, X2, .., XN, Y, PL) % % INPUTS: Y - dependent variable % X(1..N) - input variables % PL - parameter list % % OUTPUT: P - a pest object with the N+1 elements % % % PARAMETERS: % 'dy' - uncertainty on the dependent variable % 'dx' - uncertainties on the input variables % 'p0' - initial guess on the fit parameters to propagate uncertainities % in the input variables X(1..N) to the dependent variable Y % % <a href="matlab:utils.helper.displayMethodInfo('ao', 'bilinfit')">Parameters Description</a> % % VERSION: $Id: bilinfit.m,v 1.20 2011/04/08 08:56:11 hewitson Exp $ % % EXAMPLES: % % % 1) 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', 'm')); % x2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm')); % n = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm')); % c = [ao(1,plist('yunits','m/m')) ao(2,plist('yunits','m/m'))]; % y = c(1)*x1 + c(2)*x2 + n; % y.simplifyYunits; % % % Get a fit for the c coefficients and a constant term % p = bilinfit(x1, x2, y) % % % Do linear combination: using eval % pl_split = plist('times', [1 5]); % yfit = eval(p, split(x1, pl_split), split(x2, pl_split)); % % % Plot (compare data with fit) % iplot(y, yfit, plist('Linestyles', {'-','--'})) % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function varargout = bilinfit(varargin) % check if this is a call for parameters if utils.helper.isinfocall(varargin{:}) varargout{1} = getInfo(varargin{3}); return end % tell the system we are runing 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 [aos, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names); pl = utils.helper.collect_objects(varargin(:), 'plist', in_names); if nargout == 0 error('### bilinfit can not be used as a modifier method. Please give at least one output'); end if numel(aos) < 2 error('### bilinfit needs at least 2 inputs AOs'); end % combine plists pl = parse(pl, getDefaultPlist()); % extract parameters dy = find(pl, 'dy'); dx = find(pl, 'dx'); p0 = find(pl, 'p0'); % collect inputs Y = aos(end); X = aos(1:end-1); % collect inputs names argsname = aos(1).name; for jj = 2:numel(aos) argsname = [argsname ',' aos(jj).name]; end % get data from AOs x = X(:).y; y = Y.y; % vectors length N = length(y); % uncertainty on Y if isempty(dy) dy = 1; end if isa(dy, 'ao') % check units if Y.yunits ~= dy.yunits error('### Y and DY units are not compatible - %s %s', char(Y.yunits), char(dy.yunits)); end % extract values from AO dy = dy.y; end if isscalar(dy) % given a single value construct a vector dy = ones(N, 1) * dy; end % squares sigma2 = dy.^2; sigma2y_rms = sqrt(sum(sigma2)/N); % extract values for initial guess if (isa(p0, 'ao') || isa(p0, 'pest')) p0 = p0.y; end % uncertainty on X if ~isempty(dx) for k = 1:length(dx) dxi = dx(k); if ~isempty(dxi) if isa(dxi, 'ao') % check units if X(k).yunits ~= dxi.yunits error('### X and DX units are not compatible - %s %s', char(X.yunits), char(dxi.yunits)); end % extract values from AO dxi = dxi.y; end if isscalar(dxi) % given a single value construct a vector dxi = ones(N, 1) * dxi; end % squares sigma2xi = dxi.^2; % if A0 guess are not given if isempty(p0(k)) % set it to obtain equal error contribution to the Y error sigma2xi_rms = sqrt(sum(sigma2xi)/N); p0(k) = sigma2y_rms/sigma2xi_rms; end % add contribution to weights sigma2 = sigma2 + sigma2xi * p0(k)^2; end end end % constant term c = ones(N, 1); % build matrix m = [x c]; % solve [p, stdx, mse, s] = lscov(m, y, 1./sigma2); % scale errors and covariance matrix stdp = stdx ./ sqrt(mse); s = s ./ mse; % compute chi2 dof = N - length(p); chi2 = sum((y - lincom(m, p)).^2 ./ sigma2) / dof; % prepare model, units, names model = []; for kk = 1:length(p) switch kk case 1 units(kk) = simplify(Y.yunits/X(kk).yunits); model = ['P' num2str(kk) '*X' num2str(kk)]; xvar{kk} = ['X' num2str(kk)]; xunits{kk} = X(kk).yunits; case length(p) units(kk) = Y.yunits; model = [model ' + P' num2str(kk)]; otherwise units(kk) = simplify(Y.yunits/X(kk).yunits); model = [model ' + P' num2str(kk) '*X' num2str(kk)]; xvar{kk} = ['X' num2str(kk)]; xunits{kk} = X(kk).yunits; end names{kk} = ['P' num2str(kk)]; end model = smodel(plist('expression', model, ... 'params', names, ... 'values', p, ... 'xvar', xvar, ... 'xunits', xunits, ... 'yunits', Y.yunits)); % build the output pest object out = pest; out.setY(p); out.setDy(stdp); out.setCov(s); out.setChi2(chi2); out.setDof(dof); out.setNames(names{:}); out.setYunits(units); out.setModels(model); out.name = sprintf('bilinfit(%s)', argsname); out.addHistory(getInfo('None'), pl, ao_invars, [aos(:).hist]); % set procinfo object out.procinfo = plist('MSE', mse); % set outputs varargout{1} = out; end % computes linear combination function out = lincom(x, p) assert(size(x, 2) == length(p)); out = zeros(size(x, 1), 1); for k = 1:length(p) out = out + x(:,k) * p(k); end 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: bilinfit.m,v 1.20 2011/04/08 08:56:11 hewitson Exp $', sets, pl); ii.setModifier(false); ii.setArgsmin(2); end % get default plist function plout = getDefaultPlist() persistent pl; if ~exist('pl', 'var') || isempty(pl) pl = buildplist(); end plout = pl; end function pl = buildplist() % default plist for linear fitting pl = plist.MULTILINEAR_FIT_PLIST; end