Mercurial > hg > ltpda
view m-toolbox/classes/@ao/xfit.m @ 26:ce4df2e95a55 database-connection-manager
Remove LTPDARepositoryManager initialization
| author | Daniele Nicolodi <nicolodi@science.unitn.it> |
|---|---|
| date | Mon, 05 Dec 2011 16:20:06 +0100 |
| parents | f0afece42f48 |
| children |
line wrap: on
line source
% XFIT fit a function of x to data. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % DESCRIPTION: XFIT performs a non-linear fit of a function of x to data. % smodels fitting is also supported. % % ALGORITHM: XFIT does a chi-squared minimization by means of different % algorithms (see details in the default plist). Covariance matrix is also % computed from the Fisher's Information Matrix. In case the Information % Matrix is not positive-definite, uncertainties will not be stored in the % output. % % CALL: b = xfit(a, pl) % % INPUTS: a - input AO to fit to % pl - parameter list (see below) % % OUTPUTs: b - a pest object containing the best-fit parameters, % goodness-of-fit reduced chi-squared, fit degree-of-freedom % covariance matrix and uncertainties. Additional % quantities, like the Information Matrix, are contained % within the procinfo. The best-fit model can be evaluated % from pest\eval. % % <a href="matlab:utils.helper.displayMethodInfo('ao', 'xfit')">Parameters Description</a> % % EXAMPLES: % % % 1) Fit to a frequency-series % % % Create a frequency-series % datapl = plist('fsfcn', '0.01./(0.0001+f) + 5*abs(randn(size(f))) ', 'f1', 1e-5, 'f2', 5, 'nf', 1000, ... % 'xunits', 'Hz', 'yunits', 'N/Hz'); % data = ao(datapl); % data.setName; % % % Do fit % fitpl = plist('Function', 'P(1)./(P(2) + Xdata) + P(3)', ... % 'P0', [0.1 0.01 1]); % params = xfit(data, fitpl) % % % Evaluate model % BestModel = eval(params,plist('type','fsdata','xdata',data,'xfield','x')); % BestModel.setName; % % % Display results % iplot(data,BestModel) % % % 2) Fit to a noisy sine-wave % % % Create a noisy sine-wave % fs = 10; % nsecs = 500; % datapl = plist('waveform', 'Sine wave', 'f', 0.01, 'A', 0.6, 'fs', fs, 'nsecs', nsecs, ... % 'xunits', 's', 'yunits', 'm'); % sw = ao(datapl); % noise = ao(plist('tsfcn', '0.01*randn(size(t))', 'fs', fs, 'nsecs', nsecs)); % data = sw+noise; % data.setName; % % % Do fit % fitpl = plist('Function', 'P(1).*sin(2*pi*P(2).*Xdata + P(3))', ... % 'P0', [1 0.01 0]); % params = xfit(data, fitpl) % % % Evaluate model % BestModel = eval(params,plist('type','tsdata','xdata',sw,'xfield','x')); % BestModel.setName; % % % Display results % iplot(data,BestModel) % % % 3) Fit an smodel of a straight line to some data % % % Create a noisy straight-line % datapl = plist('xyfcn', '2.33 + 0.1*x + 0.01*randn(size(x))', 'x', 0:0.1:10, ... % 'xunits', 's', 'yunits', 'm'); % data = ao(datapl); % data.setName; % % % Model to fit % mdl = smodel('a + b*x'); % mdl.setXvar('x'); % mdl.setParams({'a', 'b'}, {1 2}); % % % Fit model % fitpl = plist('Function', mdl, 'P0', [1 1]); % params = xfit(data, fitpl) % % % Evaluate model % BestModel = eval(params,plist('xdata',data,'xfield','x')); % BestModel.setName; % % % Display results % iplot(data,BestModel) % % % 4) Fit a chirp-sine firstly starting from an initial guess (quite close % % to the true values) (bad convergency) and secondly by a Monte Carlo % % search (good convergency) % % % Create a noisy chirp-sine % fs = 10; % nsecs = 1000; % % % Model to fit and generate signal % mdl = smodel(plist('name', 'chirp', 'expression', 'A.*sin(2*pi*(f + f0.*t).*t + p)', ... % 'params', {'A','f','f0','p'}, 'xvar', 't', 'xunits', 's', 'yunits', 'm')); % % % signal % s = mdl.setValues({10,1e-4,1e-5,0.3}); % s.setXvals(0:1/fs:nsecs-1/fs); % signal = s.eval; % signal.setName; % % % noise % noise = ao(plist('tsfcn', '1*randn(size(t))', 'fs', fs, 'nsecs', nsecs)); % % % data % data = signal + noise; % data.setName; % % % Fit model from the starting guess % fitpl_ig = plist('Function', mdl, 'P0',[8,9e-5,9e-6,0]); % params_ig = xfit(data, fitpl_ig); % % % Evaluate model % BestModel_ig = eval(params_ig,plist('xdata',data,'xfield','x')); % BestModel_ig.setName; % % % Display results % iplot(data,BestModel_ig) % % % Fit model by a Monte Carlo search % fitpl_mc = plist('Function', mdl, ... % 'MonteCarlo', 'yes', 'Npoints', 1000, 'LB', [8,9e-5,9e-6,0], 'UB', [11,3e-4,2e-5,2*pi]); % params_mc = xfit(data, fitpl_mc) % % % Evaluate model % BestModel_mc = eval(params_mc,plist('xdata',data,'xfield','x')); % BestModel_mc.setName; % % % Display results % iplot(data,BestModel_mc) % % % 5) Multichannel fit of smodels % % % Ch.1 data % datapl = plist('xyfcn', '0.1*x + 0.01*randn(size(x))', 'x', 0:0.1:10, 'name', 'channel 1', ... % 'xunits', 'K', 'yunits', 'Pa'); % a1 = ao(datapl); % % Ch.2 data % datapl = plist('xyfcn', '2.5*x + 0.1*sin(2*pi*x) + 0.01*randn(size(x))', 'x', 0:0.1:10, 'name', 'channel 2', ... % 'xunits', 'K', 'yunits', 'T'); % a2 = ao(datapl); % % % Model to fit % mdl1 = smodel('a*x'); % mdl1.setXvar('x'); % mdl1.setParams({'a'}, {1}); % mdl1.setXunits('K'); % mdl1.setYunits('Pa'); % % mdl2 = smodel('b*x + a*sin(2*pi*x)'); % mdl2.setXvar('x'); % mdl2.setParams({'a','b'}, {1,2}); % mdl2.setXunits('K'); % mdl2.setYunits('T'); % % % Fit model % params = xfit(a1,a2, plist('Function', [mdl1,mdl2])); % % % evaluate model % b = eval(params, plist('index',1,'xdata',a1,'xfield','x')); % b.setName('fit Ch.1'); % r = a1-b; % r.setName('residuals'); % iplot(a1,b,r) % % b = eval(params, plist('index',2,'xdata',a2,'xfield','x')); % b.setName('fit Ch.2'); % r = a2-b; % r.setName('residuals'); % iplot(a2,b,r) % % <a href="matlab:utils.helper.displayMethodInfo('ao', 'xfit')">Parameters Description</a> % % VERSION: $Id: xfit.m,v 1.83 2011/05/12 07:46:29 mauro Exp $ % % HISTORY: 17-09-2009 G. Congedo % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function varargout = xfit(varargin) % global variables global Pidx Ydata weights modelFuncs dFcns hFcns lb ub Nch Ndata estimator % 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 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names); pl = utils.helper.collect_objects(varargin(:), 'plist', in_names); if nargout == 0 error('### xfit cannot be used as a modifier. Please give an output variable.'); end % combine plists pl = parse(pl, getDefaultPlist()); % Extract necessary parameters targetFcn = pl.find('Function'); P0 = pl.find('P0'); % ADDP = pl.find('ADDP'); userOpts = pl.find('OPTSET'); weights = pl.find('WEIGHTS'); FitUnc = pl.find('FitUnc'); UncMtd = pl.find('UncMtd'); linUnc = pl.find('LinUnc'); FastHess = pl.find('FastHess'); SymDiff = pl.find('SymDiff'); SymGrad = pl.find('SymGrad'); SymHess = pl.find('SymHess'); DiffOrder = pl.find('DiffOrder'); lb = pl.find('LB'); ub = pl.find('UB'); MCsearch = pl.find('MonteCarlo'); Npoints = pl.find('Npoints'); Noptims = pl.find('Noptims'); % SVDsearch = pl.find('SVD'); % nSVD = pl.find('nSVD'); Algorithm = pl.find('Algorithm'); % AdpScale = pl.find('AdpScale'); % only for function handle fitting pnames = pl.find('pnames'); % pvalues = pl.find('pvalues'); estimator = pl.find('estimator'); % Convert yes/no, true/false, etc. to booleans FitUnc = utils.prog.yes2true(FitUnc); linUnc = utils.prog.yes2true(linUnc); FastHess = utils.prog.yes2true(FastHess); SymDiff = utils.prog.yes2true(SymDiff); MCsearch = utils.prog.yes2true(MCsearch); % SVDsearch = utils.prog.yes2true(SVDsearch); % AdpScale = utils.prog.yes2true(AdpScale); % Check the fitting algorithm % AlgoCheck = ~strcmpi(Algorithm,'fminsearch') & ~strcmpi(Algorithm,'fminunc') & ~isempty(Algorithm); if isempty(Algorithm) Algorithm = 'fminsearch'; utils.helper.msg(msg.IMPORTANT, 'using %s as the fitting algorithm', upper(Algorithm)); elseif ischar(Algorithm) Algorithm = lower(Algorithm); switch Algorithm case 'fminsearch' utils.helper.msg(msg.IMPORTANT, 'using %s as the fitting algorithm', upper(Algorithm)); Algorithm = 'fminsearch'; case 'fminunc' if exist('fminunc','file')==2 utils.helper.msg(msg.IMPORTANT, 'using %s as the fitting algorithm', upper(Algorithm)); Algorithm = 'fminunc'; else error('### you must install Optimization Toolbox in order to use %s', upper(Algorithm)) end case 'fmincon' if exist('fmincon','file')==2 utils.helper.msg(msg.IMPORTANT, 'using %s as the fitting algorithm', upper(Algorithm)); Algorithm = 'fmincon'; else error('### you must install Optimization Toolbox in order to use %s', upper(Algorithm)) end case 'patternsearch' if exist('patternsearch','file')==2 utils.helper.msg(msg.IMPORTANT, 'using %s as the fitting algorithm', upper(Algorithm)); Algorithm = 'patternsearch'; else error('### you must install Genetic Algorithm and Direct Search Toolbox in order to use %s', upper(Algorithm)) end case 'ga' if exist('ga','file')==2 utils.helper.msg(msg.IMPORTANT, 'using %s as the fitting algorithm', upper(Algorithm)); Algorithm = 'ga'; else error('### you must install Genetic Algorithm and Direct Search Toolbox in order to use %s', upper(Algorithm)) end case 'simulannealbnd' if exist('simulannealbnd','file')==2 utils.helper.msg(msg.IMPORTANT, 'using %s as the fitting algorithm', upper(Algorithm)); Algorithm = 'simulannealbnd'; else error('### you must install Genetic Algorithm and Direct Search Toolbox in order to use %s', upper(Algorithm)) end otherwise error('### unknown fitting algorithm') end else error('### unknown format for ALGORITHM parameter') end % Estimator if isempty(estimator) estimator = 'chi2'; elseif ~strcmp(estimator,{'chi2','abs','log','median'}) error('### unknown name of the estimator') end % Data we will fit to Xdata = as.x; Ydata = as.y; % dYdata = as.dy; % Number of data point per each channel Ndata = numel(as(1).x); % Number of channels Nch = numel(as); multiCh = Nch-1; % Number of models if all(isa(targetFcn, 'smodel')) Nmdls = numel(targetFcn); elseif iscell(targetFcn) Nmdls = numel(targetFcn); else Nmdls = 1; end % consistency check on the data units Xunits = as.xunits; % Xunits = as(1).xunits.strs; if Nch>1 XunitsCheck = Xunits(1).strs; for kk=2:Nch if ~strcmp(Xunits(kk).strs,XunitsCheck) error('### in multi-channel fitting the xunits of all data objects must be the same') end end end Xunits = as(1).xunits; Yunits = as.yunits; % consistency check on the inputs if isempty(targetFcn) error('### please specify at least a model'); end if Nch~=Nmdls error('### number of data channels and models do not match') end for kk=1:Nch if any(numel(as(kk).x)~=Ndata) error('### the number of data is not self-consistent: please check that all channels have the same length') end end for kk=1:Nch for kkk=1:Nch if any(Xdata(:,kk)~=Xdata(:,kkk)) error('### in multi-channel fitting the x-field of all data objects must be the same') end end end cellMdl = iscell(targetFcn); if multiCh if cellMdl for kk=1:Nmdls if ~isa(targetFcn{kk}, 'function_handle') error('### please use a cell array of function handles') end end else for kk=1:Nmdls if ~isa(targetFcn(kk), 'smodel') error('### multi-channel fitting is only possible with smodels') end if isempty(targetFcn(kk).expr) error('### please specify a target function for all smodels'); end if ~isempty(targetFcn(kk).values) & size(targetFcn(kk).values)~=size(targetFcn(kk).params) error('### please give an initial value for each parameter') end end if ~isempty(P0) error('in multi-channel fitting the initial values for the parameters must be within each smodel') end checkExistAllParams = 0; for kk=1:Nch checkExistAllParams = checkExistAllParams + ~isempty(targetFcn(kk).values); end if checkExistAllParams==0 && ~MCsearch error('### please give initial values for all parameters or use Monte Carlo instead') end end end % check P0 if isempty(P0) && ~MCsearch for kk=1:Nmdls if isa(targetFcn(kk), 'smodel') && isempty(targetFcn(kk).values) error('### please give initial values for all parameters or use Monte Carlo instead') elseif ischar(targetFcn) error('### please give initial values for all parameters or use Monte Carlo instead') % elseif cellMdl % error('### please give initial values for all parameters or use Monte Carlo instead') end end end % Extract anonymous functions if multiCh if ~cellMdl % concatenate models and parameters [params,mdl_params,Pidx] = cat_mdls(targetFcn); else%if iscell(P0) [params,mdl_params,Pidx] = cat_mdls_cell(targetFcn, pnames); % else % params = pnames; % Pidx = cell(1,Nch); % mdl_params = pnames; % for ii=1:Nch % Pidx{ii} = ones(1,numel(pnames)); % end end % create the full initial value array if ~MCsearch && ~cellMdl P0 = zeros(1,numel(params)); for ii=1:numel(params) for kk=1:Nmdls for jj=1:numel(targetFcn(kk).params) if strcmp(params{ii},targetFcn(kk).params{jj}) P0(ii) = targetFcn(kk).values{jj}; end end end end elseif cellMdl if isempty(P0)% || ~iscell(pvalues) error('### please give initial values') end if isempty(pnames) || ~iscell(pnames) error('### please give parameter names in a cell array') end if size(P0)~=size(pnames) error('### the size of pnames and pvalues does not match') end % P0 = zeros(1,numel(params)); % for ii=1:numel(P0) % P0(ii) = pvalues{ii}; % end if iscell(P0) && ~MCsearch for ii=1:numel(P0) if isempty(P0{ii}) error('### please give initial values for all parameters or use Monte Carlo instead'); end end for ii=1:numel(params) for kk=1:Nmdls for jj=1:numel(pnames{kk}) if strcmp(params{ii},pnames{kk}{jj}) P0new(ii) = P0{kk}(jj); end end end end P0 = P0new; end end if all(isa(targetFcn,'smodel')) % anonymous fcns modelFuncs = cell(1,Nch); for kk=1:Nch targetFcn(kk).setXvals(Xdata(:,kk)); modelFuncs{kk} = targetFcn(kk).fitfunc; end end % if cellMdl % for kk=1:Nch % modelFuncs{kk} = targetFcn{kk}; % end % end else % Check parameters if isempty(P0) if isa(targetFcn, 'smodel') P0 = [targetFcn.values{:}]; elseif isempty(P0) && ~MCsearch error('### please give initial values for all parameters or use Monte Carlo instead'); end end if size(P0)~=[1 numel(P0)] P0 = P0'; end % convert input regular expression to anonymous function only for % single-channel fitting % create anonymouse function from user input expression if ischar(targetFcn) checkFcn = regexp(targetFcn, 'Xdata'); if isempty(checkFcn) error('### when using a string expression for the input model, the independent variable must be named as Xdata') end % tfunc = eval(['@(P,Xdata,ADDP) (',targetFcn,')']); tfunc = eval(['@(P,Xdata) (',targetFcn,')']); % now create another anonymous function that only depends on the % parameters % modelFunc = @(x)tfunc(x, Xdata, ADDP); modelFuncs{1} = @(x)tfunc(x, Xdata); elseif isa(targetFcn,'smodel') targetFcn.setXvals(Xdata); modelFuncs{1} = targetFcn.fitfunc; % elseif isa(targetFcn,'function_handle') % modelFunc = targetFcn; end end if cellMdl for kk=1:Nch modelFuncs{kk} = targetFcn{kk}; end % if ~multiCh % modelFunc = modelFuncs{1}; % end end % check lb and ub % constrained search or not conSearch = ~isempty(lb) || ~isempty(ub); if MCsearch || conSearch if isempty(lb) || isempty(ub) error('### please give LB and UB') end if multiCh && (~iscell(lb) || ~iscell(ub)) error('### in multi-channel fitting upper and lower bounds must be cell array') end if size(lb)~=size(ub) % | size(lb)~=size(P0) | size(ub)~=size(P0) error('### LB and UB must be of the same size'); end if multiCh && numel(lb)~=Nch && numel(ub)~=Nch error('### in multi-channel fitting LB and UB must be cell array whose number of elements is equal to the number of models'); end if ~multiCh && ~all(lb<=ub) error('### UB must me greater equal to LB'); end if multiCh for kk=1:Nmdls if numel(lb{kk})~=numel(mdl_params{kk}) error('### please give the proper number of values for LB for each model') end if numel(ub{kk})~=numel(mdl_params{kk}) error('### please give the proper number of values for UB for each model') end if ~all(lb{kk}<=ub{kk}) error('### UB must me greater equal to LB for all parameters and models'); end if size(lb{kk})~=[1 numel(lb{kk})] lb{kk} = lb{kk}'; ub{kk} = ub{kk}'; end end end if ~multiCh if size(lb)~=[1 numel(lb)] lb = lb'; ub = ub'; end end end % create the full bounds array if (MCsearch || conSearch) && multiCh && ~cellMdl lb_full = zeros(1,numel(params)); ub_full = zeros(1,numel(params)); for ii=1:numel(params) for kk=1:Nmdls for jj=1:numel(targetFcn(kk).params) if strcmp(params{ii},targetFcn(kk).params{jj}) lb_full(ii) = lb{kk}(jj); ub_full(ii) = ub{kk}(jj); end end end end lb = lb_full; ub = ub_full; elseif (MCsearch || conSearch) && multiCh && cellMdl lb_full = zeros(1,numel(params)); ub_full = zeros(1,numel(params)); for ii=1:numel(params) for kk=1:Nmdls for jj=1:numel(mdl_params{kk}) if strcmp(params{ii},mdl_params{kk}(jj)) lb_full(ii) = lb{kk}(jj); ub_full(ii) = ub{kk}(jj); end end end end lb = lb_full; ub = ub_full; % elseif cellMdl % lb = lb; % ub = ub; end % if ~iscell(ADDP) % ADDP = {ADDP}; % end % Get input options switch Algorithm case 'fminsearch' opts = optimset(@fminsearch); if isstruct(userOpts) opts = optimset(opts, userOpts); else for j=1:2:numel(userOpts) opts = optimset(opts, userOpts{j}, userOpts{j+1}); end end case 'fminunc' opts = optimset(@fminunc); if isstruct(userOpts) opts = optimset(opts, userOpts); else for j=1:2:numel(userOpts) opts = optimset(opts, userOpts{j}, userOpts{j+1}); end end case 'fmincon' opts = optimset(@fmincon); if isstruct(userOpts) opts = optimset(opts, userOpts); else for j=1:2:numel(userOpts) opts = optimset(opts, userOpts{j}, userOpts{j+1}); end end case 'patternsearch' opts = psoptimset(@patternsearch); if isstruct(userOpts) opts = psoptimset(opts, userOpts); else for j=1:2:numel(userOpts) opts = psoptimset(opts, userOpts{j}, userOpts{j+1}); end end case 'ga' opts = gaoptimset(@ga); if isstruct(userOpts) opts = gaoptimset(opts, userOpts); else for j=1:2:numel(userOpts) opts = gaoptimset(opts, userOpts{j}, userOpts{j+1}); end end case 'simulannealbnd' opts = saoptimset(@simulannealbnd); if isstruct(userOpts) opts = saoptimset(opts, userOpts); else for j=1:2:numel(userOpts) opts = saoptimset(opts, userOpts{j}, userOpts{j+1}); end end end % compute the right weights [weights,unweighted] = find_weights(as,weights); % define number of free parameters if ~MCsearch Nfreeparams = length(P0); else Nfreeparams = length(lb); end if Nch==1 Pidx{1}=1:Nfreeparams; % modelFuncs{1}=modelFunc; end % Check for user-supplied analytical gradient as cell-array of function % handles if ~isempty(SymGrad) dFcns = cell(Nmdls,1); for kk=1:Nmdls for ii=1:numel(SymGrad{kk}) dFcns{kk}{ii} = SymGrad{kk}{ii}; end end end % Check for user-supplied analytical hessian as cell-array of function % handles if ~isempty(SymGrad) && ~isempty(SymHess) hFcns = cell(Nmdls,1); for kk=1:Nmdls for jj=1:numel(dFcns{kk}) for ii=1:jj hFcns{kk}{ii,jj} = SymHess{kk}{ii,jj}; end end end end % If requested, compute the analytical gradient if SymDiff || linUnc && isempty(SymGrad) utils.helper.msg(msg.IMPORTANT, 'evaluating symbolic derivatives'); if ~isa(targetFcn,'smodel') error('### smodel functions must be used in order to do symbolic differentiation') end % compute symbolic 1st-order differentiation dFcnsSmodel = cell(Nmdls,1); for kk=1:Nmdls p = targetFcn(kk).params; for ii=1:numel(p) dFcnsSmodel{kk}(ii) = diff(targetFcn(kk),p{ii}); end end % extract anonymous function dFcns = cell(Nmdls,1); for kk=1:Nmdls for ii=1:numel(dFcnsSmodel{kk}) dFcns{kk}{ii} = dFcnsSmodel{kk}(ii).fitfunc; end end if DiffOrder==2; % compute symbolic 2nd-order differentiation hFcnsSmodel = cell(Nmdls,1); for kk=1:Nmdls p = targetFcn(kk).params; for jj=1:numel(p) for ii=1:jj hFcnsSmodel{kk}(ii,jj) = diff(dFcnsSmodel{kk}(ii),p{jj}); end end end % extract anonymous function hFcns = cell(Nmdls,1); for kk=1:Nmdls for jj=1:numel(dFcnsSmodel{kk}) for ii=1:jj hFcns{kk}{ii,jj} = hFcnsSmodel{kk}(ii,jj).fitfunc; end end end end end % Set optimset in order to take care of eventual symbolic differentiation if ~isempty(dFcns) && any(strcmp(Algorithm,{'fminunc','fmincon'})) opts = optimset(opts, 'GradObj', 'on'); end if ~isempty(hFcns) && any(strcmp(Algorithm,{'fminunc','fmincon'})) opts = optimset(opts, 'Hessian', 'on'); end % Start the best-fit search if MCsearch utils.helper.msg(msg.IMPORTANT, 'performing a Monte Carlo search'); % find best-fit by a Monte Carlo search [P,chi2,exitflag,output,h,MChistory] = ... find_bestfit_montecarlo(lb, ub, Npoints, Noptims, opts, Algorithm, FastHess); else utils.helper.msg(msg.IMPORTANT, 'looking for a best-fit from the initial guess'); % find best-fit starting from an initial guess [P,chi2,exitflag,output,h,chainHistory] = ... find_bestfit_guess(P0, opts, Algorithm, FastHess); end % degrees of freedom in the problem % dof = Nch * (Ndata - Nfreeparams); dof = Nch * Ndata - Nfreeparams; % redefine MChistory's column to put the reduced chi2s if MCsearch MChistory(:,1) = MChistory(:,1)./dof; end % redefine history to put the reduced chi2s if ~MCsearch chainHistory.fval = chainHistory.fval./dof; end % Confidence intervals contruction if FitUnc utils.helper.msg(msg.IMPORTANT, 'estimating confidence intervals'); % find best-fit errors [se,Sigma,Corr,I,H,errH,J,errJ] = ... find_errors(modelFuncs, P, chi2, dof, UncMtd, FastHess, h, weights, unweighted, linUnc, dFcns); % report issue on covariance matrix in case of % degeneracy/quasi-singularity/ill-conditioning, etc. posDef = all(diag(Sigma)>=0); if ~posDef % analysis of information matrix in order to cancel out un-important % parameters Inorm=I./norm(I); pnorms = zeros(1,size(Inorm,2)); for ii=1:size(I,2); pnorms(ii) = norm(Inorm(:,ii)); end [pnorms_sort,IX] = sort(pnorms,'descend'); if Nmdls>1 pnames_sort = params(IX); elseif isa(targetFcn,'smodel') && Nmdls==1%&& ~isempty(targetFcn.params) params = targetFcn.params; pnames_sort = params(IX); elseif isa(targetFcn,'function_handle') %&& ~isempty(targetFcn.params) pnames_sort = pnames(IX); end utils.helper.msg(msg.IMPORTANT, ['Information matrix is quasi-singular due to degeneracy or ill-conditioning: \n'... 'consider eliminating the parameters with low information.']); % ...'In the following, parameters are reported in descending order of information']); for ii=1:numel(pnorms_sort); if exist('pnames_sort','var') utils.helper.msg(msg.IMPORTANT, 'param %s: %g ', pnames_sort{ii}, pnorms_sort(ii)); else utils.helper.msg(msg.IMPORTANT, 'param %d: %g ', IX(ii), pnorms_sort(ii)); end end end end utils.helper.msg(msg.IMPORTANT, 'final best-fit found at reduced chi2, dof: %g %d ', chi2/dof, dof); if FitUnc && ~isempty(se) for kk=1:numel(P) utils.helper.msg(msg.IMPORTANT, 'best-fit param %d: %g +- %2.1g ', kk, P(kk), se(kk)); end else for kk=1:numel(P) utils.helper.msg(msg.IMPORTANT, 'best-fit param %d: %g ', kk, P(kk)); end end % check the existence of all variables if ~exist('se','var'); se = []; end if ~exist('Sigma','var'); Sigma = []; end if ~exist('Corr','var'); Corr = []; end if ~exist('I','var'); I = []; end if ~exist('H','var'); H = []; end if ~exist('errH','var'); errH = []; end if ~exist('J','var'); J = []; end if ~exist('errJ','var'); errJ = []; end if ~exist('Sigma','var'); Sigma = []; end if ~exist('MChistory','var'); MChistory = []; end % if ~exist('SVDhistory','var'); SVDhistory = []; end if ~exist('output','var'); output = []; end % Make output pest out = pest; if Nch>1 %exist('params') out.setNames(params); elseif Nch==1 && isa(targetFcn, 'smodel') out.setNames(targetFcn.params); end out.setY(P'); out.setDy(se'); out.setCov(Sigma); out.setCorr(Corr); out.setChi2(chi2/dof); out.setDof(dof); if ~MCsearch out.setChain([chainHistory.fval,chainHistory.x]); end % add the output best-fit models outFcn = targetFcn; if isa(targetFcn, 'smodel') && Nmdls>1 for kk=1:Nmdls p = P(Pidx{kk}); outFcn(kk).setValues(p); % outFcn(kk).setXvals(Xdata(:,1)); outFcn(kk).setXunits(Xunits); outFcn(kk).setYunits(Yunits(kk)); outFcn(kk).setName(['Best-fit model ' num2str(kk)]); end out.setModels(outFcn); elseif isa(targetFcn, 'smodel') && Nmdls==1 outFcn.setValues(P'); % outFcn.setXvals(Xdata(:,1)); outFcn.setXunits(Xunits); outFcn.setYunits(Yunits); outFcn.setName('Best-fit model'); out.setModels(outFcn); elseif ischar(targetFcn) % convert regular expression into smodel targetFcn = regexprep(targetFcn, 'Xdata', 'x'); for ii=1:Nfreeparams targetFcn = regexprep(targetFcn, ['P\(' num2str(ii) '\)'], ['P' num2str(ii)]); end outFcn = smodel((targetFcn)); pnames = cell(1,Nfreeparams); for ii=1:Nfreeparams pnames{ii} = ['P' num2str(ii)]; end outFcn.setParams(pnames, P'); outFcn.setXvar('x'); % outFcn.setXvals(Xdata); outFcn.setXunits(Xunits); outFcn.setYunits(Yunits); outFcn.setName('Best-fit model'); out.setModels(outFcn); out.setNames(pnames); end % Set Name, History and Procinfo if ~cellMdl mdlname = char(targetFcn); if numel(mdlname)>20 mdlname = mdlname(1:20); end out.name = sprintf('xfit(%s)', mdlname); else mdlname = func2str(targetFcn{1}); for kk=2:Nch mdlname = strcat(mdlname,[',' func2str(targetFcn{kk})]); end out.name = sprintf('xfit(%s)', mdlname); end out.addHistory(getInfo('None'), pl, ao_invars(:), [as(:).hist]); out.procinfo = plist('algorithm', Algorithm, 'exitflag', exitflag, 'output', output, ... 'InfoMatrix', I,... 'MChistory', MChistory); % 'SVDhistory', SVDhistory, 'res', res, 'hess', H, 'errhess', errH, 'jac', J, 'errjac', errJ); % Set outputs if nargout > 0 varargout{1} = out; end clear global Pidx Ydata weights modelFuncs dFcns hFcns lb ub scale Nch Ndata estimator end %-------------------------------------------------------------------------- % Included Functions %-------------------------------------------------------------------------- function [chi2,g,H] = fit_chi2(x) global Pidx Ydata weights modelFuncs dFcns hFcns lb ub scale estimator Nfreeparams = numel(x); Nmdls = numel(modelFuncs); Ndata = numel(Ydata(:,1)); mdldata = zeros(Ndata,Nmdls); for kk=1:Nmdls p = x(Pidx{kk}).*scale(Pidx{kk}); mdldata(:,kk) = modelFuncs{kk}(p); end res = (mdldata-Ydata).*weights; if all(x.*scale>=lb & x.*scale<=ub) % chi2 = res'*res; % chi2 = sum(diag(chi2)); if strcmp(estimator,'chi2') chi2 = sum(sum(res.^2)); elseif strcmp(estimator,'abs') chi2 = sum(sum(abs(res))); elseif strcmp(estimator,'log') chi2 = sum(sum(log(1+res.^2/2))); elseif strcmp(estimator,'median') dof = Ndata*Nmdls - Nfreeparams; res = reshape(res,numel(res),1); chi2 = median(abs(res))*dof; end else chi2 = 10e50; end if nargout > 1 % gradient required % if all(x.*scale>=lb & x.*scale<=ub) grad = cell(Nmdls,1); g = zeros(Nmdls,Nfreeparams); for kk=1:Nmdls p = x(Pidx{kk}).*scale(Pidx{kk}); Np = numel(p); grad{kk} = zeros(Ndata,Np); for ii=1:Np grad{kk}(:,ii) = dFcns{kk}{ii}(p); g(kk,Pidx{kk}(ii)) = 2.*res(:,kk)'*grad{kk}(:,ii); % dF=dFcns{kk}{ii}(p); % g(kk,Pidx{kk}(ii)) = 2.*sum(res(:,kk)'*dF); end end if Nmdls>1 g = sum(g); end % elseif any(x.*scale<lb) % g = repmat(-10e50,[1 Nfreeparams]); % elseif any(x.*scale>ub) % g = repmat(10e50,[1 Nfreeparams]); % end end if nargout > 2 % hessian required % hess = cell(Nmdls,1); H = zeros(Nmdls,Nfreeparams,Nfreeparams); for kk=1:Nmdls p = x(Pidx{kk}).*scale(Pidx{kk}); Np = numel(p); % hess{kk} = zeros(Ndata,Np); for jj=1:Np for ii=1:jj hF = hFcns{kk}{ii,jj}(p); H1 = sum(res(:,kk)'*hF); H2 = sum(weights(:,kk).*grad{kk}(:,ii).*grad{kk}(:,jj)); H(kk,Pidx{kk}(ii),Pidx{kk}(jj)) = 2*(H1+H2); end end end H = squeeze(sum(H,1)); % if Nmdls>1 % H = sum(H); % end H = H+triu(H,1)'; end end %-------------------------------------------------------------------------- % function chi2 = fit_chi2(P, ydata, weights, fcn) % % % evaluate model % mdldata = fcn(P); % % if ~isequal(size(mdldata), size(ydata)) % % ydata = ydata.'; % % end % % chi2 = sum((abs((mdldata-ydata)).^2).*weights); % % end %-------------------------------------------------------------------------- % function chi2 = fit_chi2_bounds(P, ydata, weights, fcn, LB, UB) % % % evaluate model % mdldata = fcn(P); % % if ~isequal(size(mdldata), size(ydata)) % % ydata = ydata.'; % % end % % if all(P>=LB & P<=UB) % chi2 = sum((abs((mdldata-ydata)).^2).*weights); % else % chi2 = 10e50; % end % % end %-------------------------------------------------------------------------- % function g = scaling(f, x, scale) % % g = f(x.*scale); % % end %-------------------------------------------------------------------------- function scale = find_scale(P0, LB, UB) if ~isempty(P0) Nfreeparams = numel(P0); else Nfreeparams = numel(LB); end % define parameter scale % scale = ones(1,Nfreeparams); % if ~isempty(LB) & ~isempty(UB) % for ii=1:Nfreeparams % if ~(LB(ii)==0 & UB(ii)==0) % if LB(ii)==0 | UB(ii)==0 % scale(ii) = max(abs(LB(ii)),abs(UB(ii))); % elseif abs(LB(ii))==abs(UB(ii)) % scale(ii) = abs(UB); % elseif LB(ii)~=UB(ii) % scale(ii) = sqrt(abs(LB(ii)) * abs(UB(ii))); % (abs(lb(ii)) + abs(ub(ii)))/2; % end % end % end % end if ~isempty(LB) && ~isempty(UB) scale = sqrt(abs(LB) .* abs(UB)); for ii=1:Nfreeparams if scale(ii)==0 scale(ii) = max(abs(LB(ii)),abs(UB(ii))); end end else scale = abs(P0); for ii=1:Nfreeparams if scale(ii)==0 scale(ii) = 1; end end end end %-------------------------------------------------------------------------- function [weightsOut,unweighted] = find_weights(as,weightsIn) Ydata = as.y; dYdata = as.dy; noWeights = isempty(weightsIn); noDY = isempty(dYdata); unweighted = noWeights & noDY; % check data uncertainties if any(dYdata==0) && ~noDY error('### some of the data uncertainties are zero: cannot fit') end if length(dYdata)~=length(Ydata) && ~noDY error('### length of Y fields and dY fields do not match') end if noWeights % The user did not input weights. % Looking for the dy field of the input data if noDY % Really no input from user weightsOut = ones(size(Ydata)); else % Uses uncertainties in Ydata to evaluate data weights weightsOut = 1./dYdata.^2; end % elseif isnumeric(weightsIn) % if size(weightsIn)~=size(Ydata) % weightsOut = weightsIn'; % end elseif numel(weightsIn) == 1 if isnumeric(weightsIn) weightsOut = weightsIn .* ones(size(Ydata)); elseif isa(weightsIn, 'ao') if weightsIn.yunits == as.yunits weightsOut = weightsIn.y; elseif size(weightsIn.data.getY)~=size(Ydata) error('### size of data ao and weights ao do not match') else error('### units for data and uncertainties do not match') end else error('### unknown format for weights parameter'); end else weightsOut = weightsIn; end end %-------------------------------------------------------------------------- function [params,mdl_params,paramidx] = cat_mdls(mdls) % This function concatenates the parameter names and values for % all the input models to construct the new parameter list. Nmdls = numel(mdls); % create the full parameter name array params_cat = horzcat(mdls(:).params); params_new = cell(1); for ii=1:numel(params_cat) if ~strcmp(params_cat{ii},params_new) params_new = [params_new params_cat{ii}]; end end params = params_new(2:numel(params_new)); % create the parameter list for each model for kk=1:Nmdls mdl_params{kk} = mdls(kk).params; end % create a cell array of indeces which map the parameter vector of each % model to the full one paramidx = cell(1,Nmdls); for kk=1:Nmdls for jdx=1:length(mdl_params{kk}) for idx=1:length(params) if (strcmp(params{idx}, mdl_params{kk}(jdx))) paramidx{kk}(jdx) = idx; break; else paramidx{kk}(jdx) = 0; end end end end end %-------------------------------------------------------------------------- function [params,mdl_params,paramidx] = cat_mdls_cell(mdls, pnames) % This function concatenates the parameter names and values for % all the input models to construct the new parameter list. Nmdls = numel(mdls); % create the full parameter name array params_cat = horzcat(pnames{:}); params_new = cell(1); for ii=1:numel(params_cat) if ~strcmp(params_cat{ii},params_new) params_new = [params_new params_cat{ii}]; end end params = params_new(2:numel(params_new)); % create the parameter list for each model for kk=1:Nmdls mdl_params{kk} = pnames{kk}; end % create a cell array of indeces which map the parameter vector of each % model to the full one paramidx = cell(1,Nmdls); for kk=1:Nmdls for jdx=1:length(mdl_params{kk}) for idx=1:length(params) if (strcmp(params{idx}, mdl_params{kk}(jdx))) paramidx{kk}(jdx) = idx; break; else paramidx{kk}(jdx) = 0; end end end end end %-------------------------------------------------------------------------- % function chi2 = fit_chi2_multich(P, Pidx, ydata, weights, fcns) % % Nmdls = numel(fcns); % Ndata = numel(ydata(:,1)); % % mdldata = zeros(Ndata,Nmdls); % for kk=1:Nmdls % p = P(Pidx{kk}); % mdldata(:,kk) = fcns{kk}(p); % end % % res = (ydata-mdldata).*weights; % % chi2 = res'*res; % chi2 = sum(diag(chi2)); % % end %-------------------------------------------------------------------------- % function chi2 = fit_chi2_multich_bounds(P, Pidx, ydata, weights, fcns, LB, UB) % % if all(P>=LB & P<=UB) % % Nmdls = numel(fcns); % Ndata = numel(ydata(:,1)); % % mdldata = zeros(Ndata,Nmdls); % for kk=1:Nmdls % p = P(Pidx{kk}); % mdldata(:,kk) = fcns{kk}(p); % end % % res = (ydata-mdldata).*weights; % % chi2 = res'*res; % chi2 = sum(diag(chi2)); % % else % chi2 = 10e50; % end % % end %-------------------------------------------------------------------------- function chi2_guesses = montecarlo(func, LB, UB, Npoints) Nfreeparams = numel(LB); % construct the guess matrix: each rows contain the chi2 and the % parameter guess guesses = zeros(Npoints,Nfreeparams); for jj=1:Nfreeparams guesses(:,jj) = LB(jj)+(UB(jj)-LB(jj))*rand(1,Npoints); end % evaluate the chi2 at each guess chi2_guesses = zeros(1,Nfreeparams+1); for ii=1:Npoints % check = ii-1/Npoints*100; % if rem(fix(check),10)==0 && check==fix(check) % % disp(sprintf('%d%% done',check)); % end chi2_guess = func(guesses(ii,:)); chi2_guesses = cat(1,chi2_guesses,[chi2_guess,guesses(ii,:)]); end chi2_guesses(1,:) = []; % sort the guesses for ascending chi2 chi2_guesses = sortrows(chi2_guesses); end %-------------------------------------------------------------------------- function [x,fval,exitflag,output,h,MChistory] = ... find_bestfit_montecarlo(LB, UB, Npoints, Noptims, opts, algorithm, FastHess) global scale Nfreeparams = length(LB); % check Npoints if isempty(Npoints) Npoints = 100000; end % check Noptims if isempty(Noptims) Noptims = 10; end if Npoints<Noptims error('### Npoints must be at least equal to Noptims') end % define parameter scale scale = find_scale([], LB, UB); % scale function to minimize % func = @(x)scaling(func, x, scale); % scale bounds LB = LB./scale; UB = UB./scale; % do a Monte Carlo search chi2_guesses = montecarlo(@fit_chi2, LB, UB, Npoints); % minimize over the first guesses fitresults = zeros(Noptims,Nfreeparams+1); exitflags = {}; outputs = {}; hs = {}; for ii=1:Noptims P0 = chi2_guesses(ii,2:Nfreeparams+1); if strcmpi(algorithm,'fminunc') [x,fval,exitflag,output,dummy,h] = fminunc(@fit_chi2, P0, opts); elseif strcmpi(algorithm,'fmincon') [x,fval,exitflag,output,dummy,dummy,h] = fmincon(@fit_chi2, P0, [], [], [], [], LB, UB, [], opts); elseif strcmpi(algorithm,'patternsearch') [x,fval,exitflag,output] = patternsearch(@fit_chi2, P0, [], [], [], [], LB, UB, [], opts); elseif strcmpi(algorithm,'ga') [x,fval,exitflag,output] = ga(@fit_chi2, numel(P0), [], [], [], [], LB, UB, [], opts); elseif strcmpi(algorithm,'simulannealbnd') [x,fval,exitflag,output] = simulannealbnd(@fit_chi2, P0, LB, UB, opts); else [x,fval,exitflag,output] = fminsearch(@fit_chi2, P0, opts); end fitresults(ii,1) = fval; fitresults(ii,2:Nfreeparams+1) = x; exitflags{ii} = exitflag; outputs{ii} = output; if FastHess && ~exist('h','var') h = fastHessian(@fit_chi2,x,fval); end if exist('h','var') hs{ii} = h; end end % sort the results [dummy, ix] = sort(fitresults(:,1)); fitresults = fitresults(ix,:); % % refine fit over the first bestfit % P0 = fitresults(1,2:Nfreeparams+1); % if strcmpi(algorithm,'fminunc') % [x,fval,exitflag,output] = fminunc(func, P0, opts); % elseif strcmpi(algorithm,'fmincon') % [x,fval,exitflag,output] = fmincon(func, P0, [], [], [], [], LB, UB, [], opts); % elseif strcmpi(algorithm,'patternsearch') % [x,fval,exitflag,output] = patternsearch(func, P0, [], [], [], [], LB, UB, [], opts); % elseif strcmpi(algorithm,'ga') % [x,fval,exitflag,output] = ga(func, numel(P0), [], [], [], [], LB, UB, [], opts); % ga does not improve the bestfit % elseif strcmpi(algorithm,'simulannealbnd') % [x,fval,exitflag,output] = simulannealbnd(func, P0, LB, UB, opts); % else % [x,fval,exitflag,output] = fminsearch(func, P0, opts); % end % set the best-fit fval = fitresults(1,1); x = fitresults(1,2:Nfreeparams+1); % scale best-fit x = x.*scale; % scale hessian if exist('hs','var') for kk=1:numel(hs) for ii=1:numel(x) for jj=1:numel(x) hs{kk}(ii,jj) = hs{kk}(ii,jj)/scale(ii)/scale(jj); end end end % else % h = []; end % scale fit results MChistory = fitresults; for ii=1:size(MChistory,1) MChistory(ii,2:Nfreeparams+1) = MChistory(ii,2:Nfreeparams+1).*scale; end % set the proper exitflag & output exitflags = exitflags(ix); outputs = outputs(ix); if exist('hs','var') && ~(isempty(hs)) hs = hs(ix); h = hs{1}; else h = []; end exitflag = exitflags{1}; output = outputs{1}; clear global scale end %-------------------------------------------------------------------------- function [x,fval,exitflag,output,h,outHistory] = ... find_bestfit_guess(P0, opts, algorithm, FastHess) global lb ub scale chainHistory if isempty(lb) lb = -Inf; end if isempty(ub) ub = Inf; end % define parameter scale scale = find_scale(P0, [], []); % scale initial guess P0 = P0./scale; % initialize history chainHistory.x = []; chainHistory.fval = []; opts=optimset(opts,'outputfcn',@outfun); % minimize over the initial guess if strcmpi(algorithm,'fminunc') [x,fval,exitflag,output,dummy,h] = fminunc(@fit_chi2, P0, opts); elseif strcmpi(algorithm,'fmincon') [x,fval,exitflag,output,dummy,dummy,h] = fmincon(@fit_chi2, P0, [], [], [], [], lb, ub, [], opts); elseif strcmpi(algorithm,'patternsearch') [x,fval,exitflag,output] = patternsearch(@fit_chi2, P0, [], [], [], [], lb, ub, [], opts); elseif strcmpi(algorithm,'ga') [x,fval,exitflag,output] = ga(@fit_chi2, numel(P0), [], [], [], [], lb, ub, [], opts); elseif strcmpi(algorithm,'simulannealbnd') [x,fval,exitflag,output] = simulannealbnd(@fit_chi2, P0, lb, ub, opts); else [x,fval,exitflag,output] = fminsearch(@fit_chi2, P0, opts); end if FastHess && ~exist('h','var') h = fastHessian(@fit_chi2,x,fval); end % scale best-fit x = x.*scale; % scale hessian if exist('h','var') for ii=1:numel(x) for jj=1:numel(x) h(ii,jj) = h(ii,jj)/scale(ii)/scale(jj); end end else h = []; end outHistory = chainHistory; outHistory.x = chainHistory.x.*repmat(scale,size(chainHistory.x,1),1); clear global lb ub scale chainHistory end %-------------------------------------------------------------------------- function stop = outfun(x,optimvalues,state) global chainHistory stop = false; if strcmp(state,'iter') chainHistory.fval = [chainHistory.fval; optimvalues.fval]; chainHistory.x = [chainHistory.x; x]; end end %-------------------------------------------------------------------------- function hess = fastHessian(fun,x0,fx0) % fastHessian calculates the numerical Hessian of fun evaluated at x % using finite differences. nVars = numel(x0); hess = zeros(nVars); % Define stepsize stepSize = eps^(1/4)*sign(x0).*max(abs(x0),1); % Min e Max change DiffMinChange = 1e-008; DiffMaxChange = 0.1; % Make sure step size lies within DiffMinChange and DiffMaxChange stepSize = sign(stepSize+eps).*min(max(abs(stepSize),DiffMinChange),DiffMaxChange); % Calculate the upper triangle of the finite difference Hessian element % by element, using only function values. The forward difference formula % we use is % % Hessian(i,j) = 1/(h(i)*h(j)) * [f(x+h(i)*ei+h(j)*ej) - f(x+h(i)*ei) % - f(x+h(j)*ej) + f(x)] (2) % % The 3rd term in (2) is common within each column of Hessian and thus % can be reused. We first calculate that term for each column and store % it in the row vector fplus_array. fplus_array = zeros(1,nVars); for j = 1:nVars xplus = x0; xplus(j) = x0(j) + stepSize(j); % evaluate fplus = fun(xplus); fplus_array(j) = fplus; end for i = 1:nVars % For each row, calculate the 2nd term in (4). This term is common to % the whole row and thus it can be reused within the current row: we % store it in fplus_i. xplus = x0; xplus(i) = x0(i) + stepSize(i); % evaluate fplus_i = fun(xplus); for j = i:nVars % start from i: only upper triangle % Calculate the 1st term in (2); this term is unique for each element % of Hessian and thus it cannot be reused. xplus = x0; xplus(i) = x0(i) + stepSize(i); xplus(j) = xplus(j) + stepSize(j); % evaluate fplus = fun(xplus); hess(i,j) = (fplus - fplus_i - fplus_array(j) + fx0)/(stepSize(i)*stepSize(j)); end end % Fill in the lower triangle of the Hessian hess = hess + triu(hess,1)'; end %-------------------------------------------------------------------------- % function [x,fval,exitflag,output,h,SVDhistory] = ... % find_bestfit_guess_SVD(P0, lb, ub, opts, algorithm, FastHess, nSVD) % % if isempty(lb) % lb = -Inf; % end % if isempty(ub) % ub = Inf; % end % % % define parameter scale % scale = find_scale(P0, [], []); % % % % scale function to minimize % % func = @(x)scaling(func, x, scale); % % % scale initial guess % P0 = P0./scale; % % % scale bounds % if ~isempty(lb) % lb = lb./scale; % end % if ~isempty(ub) % ub = ub./scale; % end % % % set the guess for 0-iteration % x = P0; % % Nfreeparams = numel(P0); % % fitresults = zeros(nSVD,Nfreeparams+1); % exitflags = {}; % outputs = {}; % hs = {}; % % % start SVD loop % for ii=1:nSVD % % % minimize over the old parameter space % if strcmpi(algorithm,'fminunc') % [x,fval,exitflag,output,grad,h] = fminunc(@fit_chi2, x, opts); % elseif strcmpi(algorithm,'fmincon') % [x,fval,exitflag,output,lambda,grad,h] = fmincon(@fit_chi2, x, [], [], [], [], lb, ub, [], opts); % elseif strcmpi(algorithm,'patternsearch') % [x,fval,exitflag,output] = patternsearch(@fit_chi2, P0, [], [], [], [], lb, ub, [], opts); % elseif strcmpi(algorithm,'ga') % [x,fval,exitflag,output] = ga(@fit_chi2, numel(P0), [], [], [], [], lb, ub, [], opts); % elseif strcmpi(algorithm,'simulannealbnd') % [x,fval,exitflag,output] = simulannealbnd(@fit_chi2, P0, lb, ub, opts); % else % [x,fval,exitflag,output] = fminsearch(@fit_chi2, P0, opts); % end % % if FastHess && ~exist('h','var') % h = fastHessian(func,x,fval); % end % % % compute Jacobian % J = zeros(Nch*Ndata,Np); % for kk=1:Nch % for ll=1:Np % J(((ii-1)*Ndata+1):ii*Ndata,ll) = dFcns{kk}{ll}(x.*scale); % end % end % % % take SVD decomposition of hessian % [U,S,V] = svd(J); % % % cancel out column with very low information % thresh = 100*eps; % idx = (diag(S)./norm(diag(S)))<=thresh; % pos = find(idx~=0); % % % reduced matrix % Sred = S; % Ured = U; % Vred = V; % Sred(:,pos) = []; % Sred(pos,:) = []; % Ured(:,pos) = []; % Ured(pos,:) = []; % Vred(:,pos) = []; % Vred(pos,:) = []; % % % % start from here inserting norm-rescaling... also in funcSVD! % % % guess in the new parameter space % % pay attention because here x is row-vector and xSVD ic column-vector % xSVD = U'*x'; % xSVD(pos) = []; % % % % bounds in the new parameter space % % if ~isempty(lb) && ~isempty(ub) % % lbSVD = U'*lb'; % % ubSVD = U'*ub'; % % end % % % do the change in the variables % funcSVD = @(xSVD)func_SVD(func,xSVD,U,pos); % % % minimize over the new parameter space % if strcmpi(algorithm,'fminunc') % [xSVD,fvalSVD,exitflag,output,grad,hSVD] = fminunc(funcSVD, xSVD, opts); % elseif strcmpi(algorithm,'fmincon') % [xSVD,fvalSVD,exitflag,output,lambda,grad,hSVD] = fmincon(funcSVD, xSVD, [], [], [], [], [], [], [], opts); % elseif strcmpi(algorithm,'patternsearch') % [xSVD,fvalSVD,exitflag,output] = patternsearch(funcSVD, xSVD, [], [], [], [], [], [], [], opts); % elseif strcmpi(algorithm,'ga') % [xSVD,fvalSVD,exitflag,output] = ga(funcSVD, numel(xSVD), [], [], [], [], [], [], [], opts); % elseif strcmpi(algorithm,'simulannealbnd') % [xSVD,fvalSVD,exitflag,output] = simulannealbnd(funcSVD, xSVD, [], [], opts); % else % [xSVD,fvalSVD,exitflag,output] = fminsearch(funcSVD, xSVD, opts); % end % % if FastHess && ~exist('hSVD','var') % hSVD = fastHessian(funcSVD,xSVD,fvalSVD); % end % % % if fvalSVD<fval % % % take the new % % x = (U'*xSVD)'; % % fval = fvalSVD; % % % trasformare la nuova h nella vecchia % % h = U'*hSVD; % % else % % % take the old % % end % % % return to the full parameter vector % xSVDfull = zeros(numel(xSVD)+numel(pos),1); % setVals = setdiff(1:numel(xSVDfull),pos); % xSVDfull(setVals) = xSVD; % % x = xb; % % % back to the old parameter space % x = (U*xSVDfull)'; % fval = fvalSVD; % hRed = Ured*hSVD*Vred'; % h = hRed; % for kk=1:numel(pos) % h(pos(kk),:) = zeros(1,size(h,2)); % h(:,pos(kk)) = zeros(size(h,1),1); % end % % fitresults(ii,1) = fval; % fitresults(ii,2:Nfreeparams+1) = x; % exitflags{ii} = exitflag; % outputs{ii} = output; % hs{ii} = h; % % end % % % sort the results % [chi2s, ix] = sort(fitresults(:,1)); % fitresults = fitresults(ix,:); % % % set the best-fit % fval = fitresults(1,1); % x = fitresults(1,2:Nfreeparams+1); % % % scale best-fit % x = x.*scale; % % % scale hessian % for ii=1:numel(x) % for jj=1:numel(x) % h(ii,jj) = h(ii,jj)/scale(ii)/scale(jj); % end % end % % % scale fit results % SVDhistory = fitresults; % for ii=1:size(SVDhistory,1) % SVDhistory(ii,2:Nfreeparams+1) = SVDhistory(ii,2:Nfreeparams+1).*scale; % end % % % set the proper exitflag & output % exitflags = exitflags(ix); % outputs = outputs(ix); % hs = hs(ix); % exitflag = exitflags{1}; % output = outputs{1}; % h = hs{1}; % % end %-------------------------------------------------------------------------- % function z = func_SVD(func,y,U,pos) % % yFull = zeros(numel(y)+numel(pos),1); % setVals = setdiff(1:numel(yFull),pos); % yFull(setVals) = y; % % x = (U*yFull)'; % % z = func(x); % % end %-------------------------------------------------------------------------- % function adaptive_scaling(P0, scaleIn, func, Niter, opts, algorithm) % % Nfreeparams = numel(P0); % % for ii=1:Niter % % % new guess % P0_new = P0; % % new scale % scale = find_scale(P0, [], []); % % % scale bounds % LB = LB./scale; % UB = UB./scale; % % % if ~multiCh % % rescale model % modelFunc = @(x)scaling(modelFunc, x, scale_new./scale); % % % define function to minimize % if isempty(lb) & isempty(ub) % func = @(x)fit_chi2(x, Ydata, weights, modelFunc); % else % lb = lb.*scale./scale_new; % ub = ub.*scale./scale_new; % func = @(x)fit_chi2_bounds(x, Ydata, weights, modelFunc, lb, ub); % end % % % else % % % func = @(x)fit_chi2_multich(params, x, Ydata, mdl_params, modelFuncs); % % func = @(x)scaling(func, x, scale_new/scale); % % end % % % Make new fit % [x,fval_new,exitflag,output] = fminsearch(func, P0_new./scale_new, opts); % % stopping criterion % if abs(fval_new-fval) <= 3*eps(fval_new-fval) % fval = fval_new; % scale = scale_new; % break % end % fval = fval_new; % scale = scale_new; % end % end % % end %-------------------------------------------------------------------------- function [se,Sigma,Corr,I,H,errH,J,errJ] = find_errors(modelFcn, x, fval, dof, UncMtd, FastHess, h, weights, unweighted, linUnc, dFcns) global scale Nch Ndata lb ub % set infinite bounds lb = repmat(-Inf,size(x)); ub = repmat(Inf,size(x)); % check on UncMtd and FastHess if strcmpi(UncMtd,'jacobian') FastHess = 0; elseif FastHess==1 && isempty(h) FastHess = 0; end Nfreeparams = numel(x); if ~FastHess % find new scale based on the best-fit found scale = find_scale(x, [], []); else scale = ones(1,Nfreeparams); end % scale best-fit x = x./scale; % standard deviation of the residuals sdr = sqrt(fval/dof); % compute information matrix from linear approximation: error propagation % from symbolical gradient if linUnc % compute Jacobian J = zeros(Nch*Ndata,Nfreeparams); for ii=1:Nch for ll=1:Nfreeparams J(((ii-1)*Ndata+1):ii*Ndata,ll) = dFcns{ii}{ll}(x.*scale); end end elseif ~strcmpi(UncMtd,'jacobian') || numel(modelFcn)~=1 && ~linUnc % Hessian method. % Brief description. % Here I use the most general method: Hessian of chi2 function. For further readings see % Numerical Recipes. Note that for inv, H should be non-singular. % Badly-scaled matrices may take to singularities: use Cholesky % factorization instead. % In practice, a singular Hessian matrix may go out for the following main reasons: % 1 - you have reached a false minimum (H should be strictly % positive-definite for a well-behaved minimum); % 2 - strong correlation between some parameters; % 3 - independency of the model on some parameters. % In the last two cases you should reduce the dimensionality of the % problem. % One last point discussed here in Trento is that the Hessian method reduces % to the Jacobian, when considering linear models. % Feedbacks are welcome. % Trento, Giuseppe Congedo if ~FastHess % scale function to find the hessian % func = @(x)scaling(func, x, scale); [H,errH] = hessian(@fit_chi2,x); elseif ~isempty(h) && FastHess % here the hessian is not scaled H = h; end if unweighted % information matrix based on sdr I = H / sdr^2 / 2; else % information matrix based on canonical form I = H./2; end elseif ~linUnc % % scale function to find the jacobian % modelFcn = @(x)scaling(modelFcn, x, scale); % Jacobian method [J,errJ] = jacobianest(modelFcn{1},x); end if exist('J','var') if unweighted % information matrix based on sdr I = J'*J ./ sdr^2; else % redefine the jacobian to take care of the weights if size(weights)~=size(J(:,1)) weights = weights'; end for ii=1:size(J,2) J(:,ii) = sqrt(weights).*J(:,ii); end % information matrix based on jacobian I = J'*J; end end % check positive-definiteness % [R,p] = chol(I); % posDef = p==0; % Covariance matrix from inverse Sigma = inv(I); % % Covariance matrix from pseudo-inverse % if issparse(I) % % [U,S,V] = svds(I,Nfreeparams); % % I = full(I); % I = full(I); % [U,S,V] = svd(I); % else % [U,S,V] = svd(I); % end % Sigma = V/S*U'; % Does the covariance give proper positive-definite variances? posDef = all(diag(Sigma)>=0); if posDef % Compute correlation matrix Corr = zeros(size(Sigma)); for ii=1:size(Sigma,1) for jj=1:size(Sigma,2) Corr(ii,jj) = Sigma(ii,jj) / sqrt(Sigma(ii,ii)) / sqrt(Sigma(jj,jj)); end end % Parameter standard errors se = sqrt(diag(Sigma))'; % rescale errors and covariance matrix if ~FastHess && ~linUnc % otherwise H is already scaled in find_bestfit function se = se.*scale; for ii=1:Nfreeparams for jj=1:Nfreeparams Sigma(ii,jj) = Sigma(ii,jj)*scale(ii)*scale(jj); end end end end % rescale information matrix if ~FastHess && ~linUnc % otherwise H is already scaled in find_bestfit function for ii=1:Nfreeparams for jj=1:Nfreeparams I(ii,jj) = I(ii,jj)/scale(ii)/scale(jj); end end end if ~exist('se','var'); se = []; end if ~exist('Sigma','var'); Sigma = []; end if ~exist('Corr','var'); Corr = []; end if ~exist('I','var'); I = []; end if ~exist('H','var'); H = []; end if ~exist('errH','var'); errH = []; end if ~exist('J','var'); J = []; end if ~exist('errJ','var'); errJ = []; end clear global scale 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.sigproc, '$Id: xfit.m,v 1.83 2011/05/12 07:46:29 mauro Exp $', sets, pl); ii.setModifier(false); 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(); % Function p = param({'Function', 'The function (or symbolic model <tt>smodel</tt>) of Xdata that you want to fit. For example: ''P(1)*Xdata''.'}, paramValue.EMPTY_STRING); pl.append(p); % Weights p = param({'Weights', ['An array of weights, one value per X sample.<br>'... 'By default, <tt>xfit</tt> takes data uncertainties from the dy field of the input object.'... 'If provided, <tt>weights</tt> make <tt>xfit</tt> ignore these.'... 'Otherwise <tt>weights</tt> will be considered as ones.']}, paramValue.EMPTY_STRING); pl.append(p); % P0 p = param({'P0', ['An array of starting guesses for the parameters.<br>'... 'This is not necessary if you fit with <tt>smodel</tt>s; in that case the parameter'... 'values from the <tt>smodel</tt> are taken as the initial guess.']}, paramValue.EMPTY_DOUBLE); pl.append(p); % % ADDP % p = param({'ADDP', 'A cell-array of additional parameters to pass to the function. These will not be fit.'}, {1, {'{}'}, paramValue.OPTIONAL}); % pl.append(p); % LB p = param({'LB', ['Lower bounds for the parameters.<br>'... 'This improves the convergency. Mandatory for Monte Carlo.<br>'... 'For multichannel fitting it has to be a cell-array.']}, paramValue.EMPTY_DOUBLE); pl.append(p); % UB p = param({'UB', ['Upper bounds for the parameters.<br>'... 'This improves the convergency. Mandatory for Monte Carlo.<br>'... 'For multichannel fitting it has to be a cell-array.']}, paramValue.EMPTY_DOUBLE); pl.append(p); % Algorithm p = param({'ALGORITHM', ['A string defining the fitting algorithm.<br>'... '<tt>fminunc</tt>, <tt>fmincon</tt> require ''Optimization Toolbox'' to be installed.<br>'... '<tt>patternsearch</tt>, <tt>ga</tt>, <tt>simulannealbnd</tt> require ''Genetic Algorithm and Direct Search'' to be installed.<br>']}, ... {1, {'fminsearch', 'fminunc', 'fmincon', 'patternsearch', 'ga', 'simulannealbnd'}, paramValue.SINGLE}); pl.append(p); % OPTSET p = param({'OPTSET', ['An optimisation structure to pass to the fitting algorithm.<br>'... 'See <tt>fminsearch</tt>, <tt>fminunc</tt>, <tt>fmincon</tt>, <tt>optimset</tt>, for details.<br>'... 'See <tt>patternsearch</tt>, <tt>psoptimset</tt>, for details.<br>'... 'See <tt>ga</tt>, <tt>gaoptimset</tt>, for details.<br>'... 'See <tt>simulannealbnd</tt>, <tt>saoptimset</tt>, for details.']}, paramValue.EMPTY_STRING); pl.append(p); % FitUnc p = param({'FitUnc', 'Fit parameter uncertainties or not.'}, paramValue.YES_NO); pl.append(p); % UncMtd p = param({'UncMtd', ['Choose the uncertainties estimation method.<br>'... 'For multi-channel fitting <tt>hessian</tt> is mandatory.']}, {1, {'hessian', 'jacobian'}, paramValue.SINGLE}); pl.append(p); % FastHess p = param({'FastHess', ['Choose whether or not the hessian should be estimated from a fast-forward finite differences algorithm; '... 'use this method to achieve better performance in CPU-time, but accuracy is not ensured; '... 'otherwise, the default (computer expensive, but more accurate) method will be used.']}, paramValue.NO_YES); pl.append(p); % MonteCarlo p = param({'MonteCarlo', ['Do a Monte Carlo search in the parameter space.<br>'... 'Useful when dealing with high multiplicity of local minima. May be computer-expensive.<br>'... 'Note that, if used, P0 will be ignored. It also requires to define LB and UB.']}, paramValue.NO_YES); pl.append(p); % Npoints p = param({'Npoints', 'Set the number of points in the parameter space to be extracted.'}, 100000); pl.append(p); % Noptims p = param({'Noptims', 'Set the number of optimizations to be performed after the Monte Carlo.'}, 10); pl.append(p); % estimator p = param({'estimator', ['Choose the robust local M-estimate as the statistical estimator of the fit:<br>'... ' ''chi2'' (least square) for data with normal distribution;<br>'... ' ''abs'' (mean absolute deviation) for data with double exponential distribution;<br>'... ' ''log'' (log least square) for data with Lorentzian distribution.']}, ... {1, {'chi2', 'abs', 'log'}, paramValue.SINGLE}); pl.append(p); end % END %-------------------------------------------------------------------------- % Author: John D'Errico % e-mail: woodchips@rochester.rr.com %-------------------------------------------------------------------------- % DERIVEST SUITE %-------------------------------------------------------------------------- function [der,errest,finaldelta] = derivest(fun,x0,varargin) % DERIVEST: estimate the n'th derivative of fun at x0, provide an error estimate % usage: [der,errest] = DERIVEST(fun,x0) % first derivative % usage: [der,errest] = DERIVEST(fun,x0,prop1,val1,prop2,val2,...) % % Derivest will perform numerical differentiation of an % analytical function provided in fun. It will not % differentiate a function provided as data. Use gradient % for that purpose, or differentiate a spline model. % % The methods used by DERIVEST are finite difference % approximations of various orders, coupled with a generalized % (multiple term) Romberg extrapolation. This also yields % the error estimate provided. DERIVEST uses a semi-adaptive % scheme to provide the best estimate that it can by its % automatic choice of a differencing interval. % % Finally, While I have not written this function for the % absolute maximum speed, speed was a major consideration % in the algorithmic design. Maximum accuracy was my main goal. % % % Arguments (input) % fun - function to differentiate. May be an inline function, % anonymous, or an m-file. fun will be sampled at a set % of distinct points for each element of x0. If there are % additional parameters to be passed into fun, then use of % an anonymous function is recommended. % % fun should be vectorized to allow evaluation at multiple % locations at once. This will provide the best possible % speed. IF fun is not so vectorized, then you MUST set % 'vectorized' property to 'no', so that derivest will % then call your function sequentially instead. % % Fun is assumed to return a result of the same % shape as its input x0. % % x0 - scalar, vector, or array of points at which to % differentiate fun. % % Additional inputs must be in the form of property/value pairs. % Properties are character strings. They may be shortened % to the extent that they are unambiguous. Properties are % not case sensitive. Valid property names are: % % 'DerivativeOrder', 'MethodOrder', 'Style', 'RombergTerms' % 'FixedStep', 'MaxStep' % % All properties have default values, chosen as intelligently % as I could manage. Values that are character strings may % also be unambiguously shortened. The legal values for each % property are: % % 'DerivativeOrder' - specifies the derivative order estimated. % Must be a positive integer from the set [1,2,3,4]. % % DEFAULT: 1 (first derivative of fun) % % 'MethodOrder' - specifies the order of the basic method % used for the estimation. % % For 'central' methods, must be a positive integer % from the set [2,4]. % % For 'forward' or 'backward' difference methods, % must be a positive integer from the set [1,2,3,4]. % % DEFAULT: 4 (a second order method) % % Note: higher order methods will generally be more % accurate, but may also suffere more from numerical % problems. % % Note: First order methods would usually not be % recommended. % % 'Style' - specifies the style of the basic method % used for the estimation. 'central', 'forward', % or 'backwards' difference methods are used. % % Must be one of 'Central', 'forward', 'backward'. % % DEFAULT: 'Central' % % Note: Central difference methods are usually the % most accurate, but sometiems one must not allow % evaluation in one direction or the other. % % 'RombergTerms' - Allows the user to specify the generalized % Romberg extrapolation method used, or turn it off % completely. % % Must be a positive integer from the set [0,1,2,3]. % % DEFAULT: 2 (Two Romberg terms) % % Note: 0 disables the Romberg step completely. % % 'FixedStep' - Allows the specification of a fixed step % size, preventing the adaptive logic from working. % This will be considerably faster, but not necessarily % as accurate as allowing the adaptive logic to run. % % DEFAULT: [] % % Note: If specified, 'FixedStep' will define the % maximum excursion from x0 that will be used. % % 'Vectorized' - Derivest will normally assume that your % function can be safely evaluated at multiple locations % in a single call. This would minimize the overhead of % a loop and additional function call overhead. Some % functions are not easily vectorizable, but you may % (if your matlab release is new enough) be able to use % arrayfun to accomplish the vectorization. % % When all else fails, set the 'vectorized' property % to 'no'. This will cause derivest to loop over the % successive function calls. % % DEFAULT: 'yes' % % % 'MaxStep' - Specifies the maximum excursion from x0 that % will be allowed, as a multiple of x0. % % DEFAULT: 100 % % 'StepRatio' - Derivest uses a proportionally cascaded % series of function evaluations, moving away from your % point of evaluation. The StepRatio is the ratio used % between sequential steps. % % DEFAULT: 2 % % % See the document DERIVEST.pdf for more explanation of the % algorithms behind the parameters of DERIVEST. In most cases, % I have chosen good values for these parameters, so the user % should never need to specify anything other than possibly % the DerivativeOrder. I've also tried to make my code robust % enough that it will not need much. But complete flexibility % is in there for your use. % % % Arguments: (output) % der - derivative estimate for each element of x0 % der will have the same shape as x0. % % errest - 95% uncertainty estimate of the derivative, such that % % abs(der(j) - f'(x0(j))) < erest(j) % % finaldelta - The final overall stepsize chosen by DERIVEST % % % Example usage: % First derivative of exp(x), at x == 1 % [d,e]=derivest(@(x) exp(x),1) % d = % 2.71828182845904 % % e = % 1.02015503167879e-14 % % True derivative % exp(1) % ans = % 2.71828182845905 % % Example usage: % Third derivative of x.^3+x.^4, at x = [0,1] % derivest(@(x) x.^3 + x.^4,[0 1],'deriv',3) % ans = % 6 30 % % True derivatives: [6,30] % % % See also: gradient % % % Author: John D'Errico % e-mail: woodchips@rochester.rr.com % Release: 1.0 % Release date: 12/27/2006 par.DerivativeOrder = 1; par.MethodOrder = 4; par.Style = 'central'; par.RombergTerms = 2; par.FixedStep = []; par.MaxStep = 100; par.StepRatio = 2; par.NominalStep = []; par.Vectorized = 'yes'; na = length(varargin); if (rem(na,2)==1) error 'Property/value pairs must come as PAIRS of arguments.' elseif na>0 par = parse_pv_pairs(par,varargin); end par = check_params(par); % Was fun a string, or an inline/anonymous function? if (nargin<1) help derivest return elseif isempty(fun) error 'fun was not supplied.' elseif ischar(fun) % a character function name fun = str2func(fun); end % no default for x0 if (nargin<2) || isempty(x0) error 'x0 was not supplied' end par.NominalStep = max(x0,0.02); % was a single point supplied? nx0 = size(x0); n = prod(nx0); % Set the steps to use. if isempty(par.FixedStep) % Basic sequence of steps, relative to a stepsize of 1. delta = par.MaxStep*par.StepRatio .^(0:-1:-25)'; ndel = length(delta); else % Fixed, user supplied absolute sequence of steps. ndel = 3 + ceil(par.DerivativeOrder/2) + ... par.MethodOrder + par.RombergTerms; if par.Style(1) == 'c' ndel = ndel - 2; end delta = par.FixedStep*par.StepRatio .^(-(0:(ndel-1)))'; end % generate finite differencing rule in advance. % The rule is for a nominal unit step size, and will % be scaled later to reflect the local step size. fdarule = 1; switch par.Style case 'central' % for central rules, we will reduce the load by an % even or odd transformation as appropriate. if par.MethodOrder==2 switch par.DerivativeOrder case 1 % the odd transformation did all the work fdarule = 1; case 2 % the even transformation did all the work fdarule = 2; case 3 % the odd transformation did most of the work, but % we need to kill off the linear term fdarule = [0 1]/fdamat(par.StepRatio,1,2); case 4 % the even transformation did most of the work, but % we need to kill off the quadratic term fdarule = [0 1]/fdamat(par.StepRatio,2,2); end else % a 4th order method. We've already ruled out the 1st % order methods since these are central rules. switch par.DerivativeOrder case 1 % the odd transformation did most of the work, but % we need to kill off the cubic term fdarule = [1 0]/fdamat(par.StepRatio,1,2); case 2 % the even transformation did most of the work, but % we need to kill off the quartic term fdarule = [1 0]/fdamat(par.StepRatio,2,2); case 3 % the odd transformation did much of the work, but % we need to kill off the linear & quintic terms fdarule = [0 1 0]/fdamat(par.StepRatio,1,3); case 4 % the even transformation did much of the work, but % we need to kill off the quadratic and 6th order terms fdarule = [0 1 0]/fdamat(par.StepRatio,2,3); end end case {'forward' 'backward'} % These two cases are identical, except at the very end, % where a sign will be introduced. % No odd/even trans, but we already dropped % off the constant term if par.MethodOrder==1 if par.DerivativeOrder==1 % an easy one fdarule = 1; else % 2:4 v = zeros(1,par.DerivativeOrder); v(par.DerivativeOrder) = 1; fdarule = v/fdamat(par.StepRatio,0,par.DerivativeOrder); end else % par.MethodOrder methods drop off the lower order terms, % plus terms directly above DerivativeOrder v = zeros(1,par.DerivativeOrder + par.MethodOrder - 1); v(par.DerivativeOrder) = 1; fdarule = v/fdamat(par.StepRatio,0,par.DerivativeOrder+par.MethodOrder-1); end % correct sign for the 'backward' rule if par.Style(1) == 'b' fdarule = -fdarule; end end % switch on par.style (generating fdarule) nfda = length(fdarule); % will we need fun(x0)? if (rem(par.DerivativeOrder,2) == 0) || ~strncmpi(par.Style,'central',7) if strcmpi(par.Vectorized,'yes') f_x0 = fun(x0); else % not vectorized, so loop f_x0 = zeros(size(x0)); for j = 1:numel(x0) f_x0(j) = fun(x0(j)); end end else f_x0 = []; end % Loop over the elements of x0, reducing it to % a scalar problem. Sorry, vectorization is not % complete here, but this IS only a single loop. der = zeros(nx0); errest = der; finaldelta = der; for i = 1:n x0i = x0(i); h = par.NominalStep(i); % a central, forward or backwards differencing rule? % f_del is the set of all the function evaluations we % will generate. For a central rule, it will have the % even or odd transformation built in. if par.Style(1) == 'c' % A central rule, so we will need to evaluate % symmetrically around x0i. if strcmpi(par.Vectorized,'yes') f_plusdel = fun(x0i+h*delta); f_minusdel = fun(x0i-h*delta); else % not vectorized, so loop f_minusdel = zeros(size(delta)); f_plusdel = zeros(size(delta)); for j = 1:numel(delta) f_plusdel(j) = fun(x0i+h*delta(j)); f_minusdel(j) = fun(x0i-h*delta(j)); end end if ismember(par.DerivativeOrder,[1 3]) % odd transformation f_del = (f_plusdel - f_minusdel)/2; else f_del = (f_plusdel + f_minusdel)/2 - f_x0(i); end elseif par.Style(1) == 'f' % forward rule % drop off the constant only if strcmpi(par.Vectorized,'yes') f_del = fun(x0i+h*delta) - f_x0(i); else % not vectorized, so loop f_del = zeros(size(delta)); for j = 1:numel(delta) f_del(j) = fun(x0i+h*delta(j)) - f_x0(i); end end else % backward rule % drop off the constant only if strcmpi(par.Vectorized,'yes') f_del = fun(x0i-h*delta) - f_x0(i); else % not vectorized, so loop f_del = zeros(size(delta)); for j = 1:numel(delta) f_del(j) = fun(x0i-h*delta(j)) - f_x0(i); end end end % check the size of f_del to ensure it was properly vectorized. f_del = f_del(:); if length(f_del)~=ndel error 'fun did not return the correct size result (fun must be vectorized)' end % Apply the finite difference rule at each delta, scaling % as appropriate for delta and the requested DerivativeOrder. % First, decide how many of these estimates we will end up with. ne = ndel + 1 - nfda - par.RombergTerms; % Form the initial derivative estimates from the chosen % finite difference method. der_init = vec2mat(f_del,ne,nfda)*fdarule.'; % scale to reflect the local delta der_init = der_init(:)./(h*delta(1:ne)).^par.DerivativeOrder; % Each approximation that results is an approximation % of order par.DerivativeOrder to the desired derivative. % Additional (higher order, even or odd) terms in the % Taylor series also remain. Use a generalized (multi-term) % Romberg extrapolation to improve these estimates. switch par.Style case 'central' rombexpon = 2*(1:par.RombergTerms) + par.MethodOrder - 2; otherwise rombexpon = (1:par.RombergTerms) + par.MethodOrder - 1; end [der_romb,errors] = rombextrap(par.StepRatio,der_init,rombexpon); % Choose which result to return % first, trim off the if isempty(par.FixedStep) % trim off the estimates at each end of the scale nest = length(der_romb); switch par.DerivativeOrder case {1 2} trim = [1 2 nest-1 nest]; case 3 trim = [1:4 nest+(-3:0)]; case 4 trim = [1:6 nest+(-5:0)]; end [der_romb,tags] = sort(der_romb); der_romb(trim) = []; tags(trim) = []; errors = errors(tags); trimdelta = delta(tags); [errest(i),ind] = min(errors); finaldelta(i) = h*trimdelta(ind); der(i) = der_romb(ind); else [errest(i),ind] = min(errors); finaldelta(i) = h*delta(ind); der(i) = der_romb(ind); end end end % mainline end function [jac,err] = jacobianest(fun,x0) % gradest: estimate of the Jacobian matrix of a vector valued function of n variables % usage: [jac,err] = jacobianest(fun,x0) % % % arguments: (input) % fun - (vector valued) analytical function to differentiate. % fun must be a function of the vector or array x0. % % x0 - vector location at which to differentiate fun % If x0 is an nxm array, then fun is assumed to be % a function of n*m variables. % % % arguments: (output) % jac - array of first partial derivatives of fun. % Assuming that x0 is a vector of length p % and fun returns a vector of length n, then % jac will be an array of size (n,p) % % err - vector of error estimates corresponding to % each partial derivative in jac. % % % Example: (nonlinear least squares) % xdata = (0:.1:1)'; % ydata = 1+2*exp(0.75*xdata); % fun = @(c) ((c(1)+c(2)*exp(c(3)*xdata)) - ydata).^2; % % [jac,err] = jacobianest(fun,[1 1 1]) % % jac = % -2 -2 0 % -2.1012 -2.3222 -0.23222 % -2.2045 -2.6926 -0.53852 % -2.3096 -3.1176 -0.93528 % -2.4158 -3.6039 -1.4416 % -2.5225 -4.1589 -2.0795 % -2.629 -4.7904 -2.8742 % -2.7343 -5.5063 -3.8544 % -2.8374 -6.3147 -5.0518 % -2.9369 -7.2237 -6.5013 % -3.0314 -8.2403 -8.2403 % % err = % 5.0134e-15 5.0134e-15 0 % 5.0134e-15 0 2.8211e-14 % 5.0134e-15 8.6834e-15 1.5804e-14 % 0 7.09e-15 3.8227e-13 % 5.0134e-15 5.0134e-15 7.5201e-15 % 5.0134e-15 1.0027e-14 2.9233e-14 % 5.0134e-15 0 6.0585e-13 % 5.0134e-15 1.0027e-14 7.2673e-13 % 5.0134e-15 1.0027e-14 3.0495e-13 % 5.0134e-15 1.0027e-14 3.1707e-14 % 5.0134e-15 2.0053e-14 1.4013e-12 % % (At [1 2 0.75], jac should be numerically zero) % % % See also: derivest, gradient, gradest % % % Author: John D'Errico % e-mail: woodchips@rochester.rr.com % Release: 1.0 % Release date: 3/6/2007 % get the length of x0 for the size of jac nx = numel(x0); MaxStep = 100; StepRatio = 2; % was a string supplied? if ischar(fun) fun = str2func(fun); end % get fun at the center point f0 = fun(x0); f0 = f0(:); n = length(f0); if n==0 % empty begets empty jac = zeros(0,nx); err = jac; return end relativedelta = MaxStep*StepRatio .^(0:-1:-25); nsteps = length(relativedelta); % total number of derivatives we will need to take jac = zeros(n,nx); err = jac; for i = 1:nx x0_i = x0(i); if x0_i ~= 0 delta = x0_i*relativedelta; else delta = relativedelta; end % evaluate at each step, centered around x0_i % difference to give a second order estimate fdel = zeros(n,nsteps); for j = 1:nsteps fdif = fun(swapelement(x0,i,x0_i + delta(j))) - ... fun(swapelement(x0,i,x0_i - delta(j))); fdel(:,j) = fdif(:); end % these are pure second order estimates of the % first derivative, for each trial delta. derest = fdel.*repmat(0.5 ./ delta,n,1); % The error term on these estimates has a second order % component, but also some 4th and 6th order terms in it. % Use Romberg exrapolation to improve the estimates to % 6th order, as well as to provide the error estimate. % loop here, as rombextrap coupled with the trimming % will get complicated otherwise. for j = 1:n [der_romb,errest] = rombextrap(StepRatio,derest(j,:),[2 4]); % trim off 3 estimates at each end of the scale nest = length(der_romb); trim = [1:3, nest+(-2:0)]; [der_romb,tags] = sort(der_romb); der_romb(trim) = []; tags(trim) = []; errest = errest(tags); % now pick the estimate with the lowest predicted error [err(j,i),ind] = min(errest); jac(j,i) = der_romb(ind); end end end % mainline function end function [grad,err,finaldelta] = gradest(fun,x0) % gradest: estimate of the gradient vector of an analytical function of n variables % usage: [grad,err,finaldelta] = gradest(fun,x0) % % Uses derivest to provide both derivative estimates % and error estimates. fun needs not be vectorized. % % arguments: (input) % fun - analytical function to differentiate. fun must % be a function of the vector or array x0. % % x0 - vector location at which to differentiate fun % If x0 is an nxm array, then fun is assumed to be % a function of n*m variables. % % arguments: (output) % grad - vector of first partial derivatives of fun. % grad will be a row vector of length numel(x0). % % err - vector of error estimates corresponding to % each partial derivative in grad. % % finaldelta - vector of final step sizes chosen for % each partial derivative. % % % Example: % [grad,err] = gradest(@(x) sum(x.^2),[1 2 3]) % grad = % 2 4 6 % err = % 5.8899e-15 1.178e-14 0 % % % Example: % At [x,y] = [1,1], compute the numerical gradient % of the function sin(x-y) + y*exp(x) % % z = @(xy) sin(diff(xy)) + xy(2)*exp(xy(1)) % % [grad,err ] = gradest(z,[1 1]) % grad = % 1.7183 3.7183 % err = % 7.537e-14 1.1846e-13 % % % Example: % At the global minimizer (1,1) of the Rosenbrock function, % compute the gradient. It should be essentially zero. % % rosen = @(x) (1-x(1)).^2 + 105*(x(2)-x(1).^2).^2; % [g,err] = gradest(rosen,[1 1]) % g = % 1.0843e-20 0 % err = % 1.9075e-18 0 % % % See also: derivest, gradient % % % Author: John D'Errico % e-mail: woodchips@rochester.rr.com % Release: 1.0 % Release date: 2/9/2007 % get the size of x0 so we can reshape % later. sx = size(x0); % total number of derivatives we will need to take nx = numel(x0); grad = zeros(1,nx); err = grad; finaldelta = grad; for ind = 1:nx [grad(ind),err(ind),finaldelta(ind)] = derivest( ... @(xi) fun(swapelement(x0,ind,xi)), ... x0(ind),'deriv',1,'vectorized','no', ... 'methodorder',2); end end % mainline function end function [HD,err,finaldelta] = hessdiag(fun,x0) % HESSDIAG: diagonal elements of the Hessian matrix (vector of second partials) % usage: [HD,err,finaldelta] = hessdiag(fun,x0) % % When all that you want are the diagonal elements of the hessian % matrix, it will be more efficient to call HESSDIAG than HESSIAN. % HESSDIAG uses DERIVEST to provide both second derivative estimates % and error estimates. fun needs not be vectorized. % % arguments: (input) % fun - SCALAR analytical function to differentiate. % fun must be a function of the vector or array x0. % % x0 - vector location at which to differentiate fun % If x0 is an nxm array, then fun is assumed to be % a function of n*m variables. % % arguments: (output) % HD - vector of second partial derivatives of fun. % These are the diagonal elements of the Hessian % matrix, evaluated at x0. % HD will be a row vector of length numel(x0). % % err - vector of error estimates corresponding to % each second partial derivative in HD. % % finaldelta - vector of final step sizes chosen for % each second partial derivative. % % % Example usage: % [HD,err] = hessdiag(@(x) x(1) + x(2)^2 + x(3)^3,[1 2 3]) % HD = % 0 2 18 % % err = % 0 0 0 % % % See also: derivest, gradient, gradest % % % Author: John D'Errico % e-mail: woodchips@rochester.rr.com % Release: 1.0 % Release date: 2/9/2007 % get the size of x0 so we can reshape % later. sx = size(x0); % total number of derivatives we will need to take nx = numel(x0); HD = zeros(1,nx); err = HD; finaldelta = HD; for ind = 1:nx [HD(ind),err(ind),finaldelta(ind)] = derivest( ... @(xi) fun(swapelement(x0,ind,xi)), ... x0(ind),'deriv',2,'vectorized','no'); end end % mainline function end function [hess,err] = hessian(fun,x0) % hessian: estimate elements of the Hessian matrix (array of 2nd partials) % usage: [hess,err] = hessian(fun,x0) % % Hessian is NOT a tool for frequent use on an expensive % to evaluate objective function, especially in a large % number of dimensions. Its computation will use roughly % O(6*n^2) function evaluations for n parameters. % % arguments: (input) % fun - SCALAR analytical function to differentiate. % fun must be a function of the vector or array x0. % fun does not need to be vectorized. % % x0 - vector location at which to compute the Hessian. % % arguments: (output) % hess - nxn symmetric array of second partial derivatives % of fun, evaluated at x0. % % err - nxn array of error estimates corresponding to % each second partial derivative in hess. % % % Example usage: % Rosenbrock function, minimized at [1,1] % rosen = @(x) (1-x(1)).^2 + 105*(x(2)-x(1).^2).^2; % % [h,err] = hessian(rosen,[1 1]) % h = % 842 -420 % -420 210 % err = % 1.0662e-12 4.0061e-10 % 4.0061e-10 2.6654e-13 % % % Example usage: % cos(x-y), at (0,0) % Note: this hessian matrix will be positive semi-definite % % hessian(@(xy) cos(xy(1)-xy(2)),[0 0]) % ans = % -1 1 % 1 -1 % % % See also: derivest, gradient, gradest, hessdiag % % % Author: John D'Errico % e-mail: woodchips@rochester.rr.com % Release: 1.0 % Release date: 2/10/2007 % parameters that we might allow to change params.StepRatio = 2; params.RombergTerms = 3; % get the size of x0 so we can reshape % later. sx = size(x0); % was a string supplied? if ischar(fun) fun = str2func(fun); end % total number of derivatives we will need to take nx = length(x0); % get the diagonal elements of the hessian (2nd partial % derivatives wrt each variable.) [hess,err] = hessdiag(fun,x0); % form the eventual hessian matrix, stuffing only % the diagonals for now. hess = diag(hess); err = diag(err); if nx<2 % the hessian matrix is 1x1. all done return end % get the gradient vector. This is done only to decide % on intelligent step sizes for the mixed partials [grad,graderr,stepsize] = gradest(fun,x0); % Get params.RombergTerms+1 estimates of the upper % triangle of the hessian matrix dfac = params.StepRatio.^(-(0:params.RombergTerms)'); for i = 2:nx for j = 1:(i-1) dij = zeros(params.RombergTerms+1,1); for k = 1:(params.RombergTerms+1) dij(k) = fun(x0 + swap2(zeros(sx),i, ... dfac(k)*stepsize(i),j,dfac(k)*stepsize(j))) + ... fun(x0 + swap2(zeros(sx),i, ... -dfac(k)*stepsize(i),j,-dfac(k)*stepsize(j))) - ... fun(x0 + swap2(zeros(sx),i, ... dfac(k)*stepsize(i),j,-dfac(k)*stepsize(j))) - ... fun(x0 + swap2(zeros(sx),i, ... -dfac(k)*stepsize(i),j,dfac(k)*stepsize(j))); end dij = dij/4/prod(stepsize([i,j])); dij = dij./(dfac.^2); % Romberg extrapolation step [hess(i,j),err(i,j)] = rombextrap(params.StepRatio,dij,[2 4]); hess(j,i) = hess(i,j); err(j,i) = err(i,j); end end end % mainline function end % ============================================ % subfunction - romberg extrapolation % ============================================ function [der_romb,errest] = rombextrap(StepRatio,der_init,rombexpon) % do romberg extrapolation for each estimate % % StepRatio - Ratio decrease in step % der_init - initial derivative estimates % rombexpon - higher order terms to cancel using the romberg step % % der_romb - derivative estimates returned % errest - error estimates % amp - noise amplification factor due to the romberg step srinv = 1/StepRatio; % do nothing if no romberg terms nexpon = length(rombexpon); rmat = ones(nexpon+2,nexpon+1); switch nexpon case 0 % rmat is simple: ones(2,1) case 1 % only one romberg term rmat(2,2) = srinv^rombexpon; rmat(3,2) = srinv^(2*rombexpon); case 2 % two romberg terms rmat(2,2:3) = srinv.^rombexpon; rmat(3,2:3) = srinv.^(2*rombexpon); rmat(4,2:3) = srinv.^(3*rombexpon); case 3 % three romberg terms rmat(2,2:4) = srinv.^rombexpon; rmat(3,2:4) = srinv.^(2*rombexpon); rmat(4,2:4) = srinv.^(3*rombexpon); rmat(5,2:4) = srinv.^(4*rombexpon); end % qr factorization used for the extrapolation as well % as the uncertainty estimates [qromb,rromb] = qr(rmat,0); % the noise amplification is further amplified by the Romberg step. % amp = cond(rromb); % this does the extrapolation to a zero step size. ne = length(der_init); rhs = vec2mat(der_init,nexpon+2,max(1,ne - (nexpon+2))); rombcoefs = rromb\(qromb.'*rhs); der_romb = rombcoefs(1,:).'; % uncertainty estimate of derivative prediction s = sqrt(sum((rhs - rmat*rombcoefs).^2,1)); rinv = rromb\eye(nexpon+1); cov1 = sum(rinv.^2,2); % 1 spare dof errest = s.'*12.7062047361747*sqrt(cov1(1)); end % rombextrap % ============================================ % subfunction - vec2mat % ============================================ function mat = vec2mat(vec,n,m) % forms the matrix M, such that M(i,j) = vec(i+j-1) [i,j] = ndgrid(1:n,0:m-1); ind = i+j; mat = vec(ind); if n==1 mat = mat.'; end end % vec2mat % ============================================ % subfunction - fdamat % ============================================ function mat = fdamat(sr,parity,nterms) % Compute matrix for fda derivation. % parity can be % 0 (one sided, all terms included but zeroth order) % 1 (only odd terms included) % 2 (only even terms included) % nterms - number of terms % sr is the ratio between successive steps srinv = 1./sr; switch parity case 0 % single sided rule [i,j] = ndgrid(1:nterms); c = 1./factorial(1:nterms); mat = c(j).*srinv.^((i-1).*j); case 1 % odd order derivative [i,j] = ndgrid(1:nterms); c = 1./factorial(1:2:(2*nterms)); mat = c(j).*srinv.^((i-1).*(2*j-1)); case 2 % even order derivative [i,j] = ndgrid(1:nterms); c = 1./factorial(2:2:(2*nterms)); mat = c(j).*srinv.^((i-1).*(2*j)); end end % fdamat % ============================================ % subfunction - check_params % ============================================ function par = check_params(par) % check the parameters for acceptability % % Defaults % par.DerivativeOrder = 1; % par.MethodOrder = 2; % par.Style = 'central'; % par.RombergTerms = 2; % par.FixedStep = []; % DerivativeOrder == 1 by default if isempty(par.DerivativeOrder) par.DerivativeOrder = 1; else if (length(par.DerivativeOrder)>1) || ~ismember(par.DerivativeOrder,1:4) error 'DerivativeOrder must be scalar, one of [1 2 3 4].' end end % MethodOrder == 2 by default if isempty(par.MethodOrder) par.MethodOrder = 2; else if (length(par.MethodOrder)>1) || ~ismember(par.MethodOrder,[1 2 3 4]) error 'MethodOrder must be scalar, one of [1 2 3 4].' elseif ismember(par.MethodOrder,[1 3]) && (par.Style(1)=='c') error 'MethodOrder==1 or 3 is not possible with central difference methods' end end % style is char valid = {'central', 'forward', 'backward'}; if isempty(par.Style) par.Style = 'central'; elseif ~ischar(par.Style) error 'Invalid Style: Must be character' end ind = find(strncmpi(par.Style,valid,length(par.Style))); if (length(ind)==1) par.Style = valid{ind}; else error(['Invalid Style: ',par.Style]) end % vectorized is char valid = {'yes', 'no'}; if isempty(par.Vectorized) par.Vectorized = 'yes'; elseif ~ischar(par.Vectorized) error 'Invalid Vectorized: Must be character' end ind = find(strncmpi(par.Vectorized,valid,length(par.Vectorized))); if (length(ind)==1) par.Vectorized = valid{ind}; else error(['Invalid Vectorized: ',par.Vectorized]) end % RombergTerms == 2 by default if isempty(par.RombergTerms) par.RombergTerms = 2; else if (length(par.RombergTerms)>1) || ~ismember(par.RombergTerms,0:3) error 'RombergTerms must be scalar, one of [0 1 2 3].' end end % FixedStep == [] by default if (length(par.FixedStep)>1) || (~isempty(par.FixedStep) && (par.FixedStep<=0)) error 'FixedStep must be empty or a scalar, >0.' end % MaxStep == 10 by default if isempty(par.MaxStep) par.MaxStep = 10; elseif (length(par.MaxStep)>1) || (par.MaxStep<=0) error 'MaxStep must be empty or a scalar, >0.' end end % check_params % ============================================ % Included subfunction - parse_pv_pairs % ============================================ function params=parse_pv_pairs(params,pv_pairs) % parse_pv_pairs: parses sets of property value pairs, allows defaults % usage: params=parse_pv_pairs(default_params,pv_pairs) % % arguments: (input) % default_params - structure, with one field for every potential % property/value pair. Each field will contain the default % value for that property. If no default is supplied for a % given property, then that field must be empty. % % pv_array - cell array of property/value pairs. % Case is ignored when comparing properties to the list % of field names. Also, any unambiguous shortening of a % field/property name is allowed. % % arguments: (output) % params - parameter struct that reflects any updated property/value % pairs in the pv_array. % % Example usage: % First, set default values for the parameters. Assume we % have four parameters that we wish to use optionally in % the function examplefun. % % - 'viscosity', which will have a default value of 1 % - 'volume', which will default to 1 % - 'pie' - which will have default value 3.141592653589793 % - 'description' - a text field, left empty by default % % The first argument to examplefun is one which will always be % supplied. % % function examplefun(dummyarg1,varargin) % params.Viscosity = 1; % params.Volume = 1; % params.Pie = 3.141592653589793 % % params.Description = ''; % params=parse_pv_pairs(params,varargin); % params % % Use examplefun, overriding the defaults for 'pie', 'viscosity' % and 'description'. The 'volume' parameter is left at its default. % % examplefun(rand(10),'vis',10,'pie',3,'Description','Hello world') % % params = % Viscosity: 10 % Volume: 1 % Pie: 3 % Description: 'Hello world' % % Note that capitalization was ignored, and the property 'viscosity' % was truncated as supplied. Also note that the order the pairs were % supplied was arbitrary. npv = length(pv_pairs); n = npv/2; if n~=floor(n) error 'Property/value pairs must come in PAIRS.' end if n<=0 % just return the defaults return end if ~isstruct(params) error 'No structure for defaults was supplied' end % there was at least one pv pair. process any supplied propnames = fieldnames(params); lpropnames = lower(propnames); for i=1:n p_i = lower(pv_pairs{2*i-1}); v_i = pv_pairs{2*i}; ind = strmatch(p_i,lpropnames,'exact'); if isempty(ind) ind = find(strncmp(p_i,lpropnames,length(p_i))); if isempty(ind) error(['No matching property found for: ',pv_pairs{2*i-1}]) elseif length(ind)>1 error(['Ambiguous property name: ',pv_pairs{2*i-1}]) end end p_i = propnames{ind}; % override the corresponding default in params params = setfield(params,p_i,v_i); %#ok end end % parse_pv_pairs % ======================================= % sub-functions % ======================================= function vec = swapelement(vec,ind,val) % swaps val as element ind, into the vector vec vec(ind) = val; end % sub-function end % ======================================= % sub-functions % ======================================= function vec = swap2(vec,ind1,val1,ind2,val2) % swaps val as element ind, into the vector vec vec(ind1) = val1; vec(ind2) = val2; end % sub-function end % End of DERIVEST SUITE