--- /dev/null
+%% Copyright (C) 1992-1994 Richard Shrager
+%% Copyright (C) 1992-1994 Arthur Jutan <jutan@charon.engga.uwo.ca>
+%% Copyright (C) 1992-1994 Ray Muzic <rfm2@ds2.uh.cwru.edu>
+%% Copyright (C) 2010, 2011 Olaf Till <i7tiol@t-online.de>
+%%
+%% This program is free software; you can redistribute it and/or modify it under
+%% the terms of the GNU General Public License as published by the Free Software
+%% Foundation; either version 3 of the License, or (at your option) any later
+%% version.
+%%
+%% This program is distributed in the hope that it will be useful, but WITHOUT
+%% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+%% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+%% details.
+%%
+%% You should have received a copy of the GNU General Public License along with
+%% this program; if not, see <http://www.gnu.org/licenses/>.
+
+%%function [f,p,cvg,iter,corp,covp,covr,stdresid,Z,r2]=
+%% leasqr(x,y,pin,F,{stol,niter,wt,dp,dFdp,options})
+%%
+%% Levenberg-Marquardt nonlinear regression of f(x,p) to y(x).
+%%
+%% Version 3.beta
+%% Optional parameters are in braces {}.
+%% x = vector or matrix of independent variables.
+%% y = vector or matrix of observed values.
+%% wt = statistical weights (same dimensions as y). These should be
+%% set to be proportional to (sqrt of var(y))^-1; (That is, the
+%% covariance matrix of the data is assumed to be proportional to
+%% diagonal with diagonal equal to (wt.^2)^-1. The constant of
+%% proportionality will be estimated.); default = ones( size (y)).
+%% pin = vec of initial parameters to be adjusted by leasqr.
+%% dp = fractional increment of p for numerical partial derivatives;
+%% default = .001*ones(size(pin))
+%% dp(j) > 0 means central differences on j-th parameter p(j).
+%% dp(j) < 0 means one-sided differences on j-th parameter p(j).
+%% dp(j) = 0 holds p(j) fixed i.e. leasqr wont change initial guess: pin(j)
+%% F = name of function in quotes or function handle; the function
+%% shall be of the form y=f(x,p), with y, x, p of the form y, x, pin
+%% as described above.
+%% dFdp = name of partial derivative function in quotes or function
+%% handle; default is 'dfdp', a slow but general partial derivatives
+%% function; the function shall be of the form
+%% prt=dfdp(x,f,p,dp,F[,bounds]). For backwards compatibility, the
+%% function will only be called with an extra 'bounds' argument if the
+%% 'bounds' option is explicitely specified to leasqr (see dfdp.m).
+%% stol = scalar tolerance on fractional improvement in scalar sum of
+%% squares = sum((wt.*(y-f))^2); default stol = .0001;
+%% niter = scalar maximum number of iterations; default = 20;
+%% options = structure, currently recognized fields are 'fract_prec',
+%% 'max_fract_change', 'inequc', 'bounds', and 'equc'. For backwards
+%% compatibility, 'options' can also be a matrix whose first and
+%% second column contains the values of 'fract_prec' and
+%% 'max_fract_change', respectively.
+%% Field 'options.fract_prec': column vector (same length as 'pin')
+%% of desired fractional precisions in parameter estimates.
+%% Iterations are terminated if change in parameter vector (chg)
+%% relative to current parameter estimate is less than their
+%% corresponding elements in 'options.fract_prec' [ie. all (abs
+%% (chg) < abs (options.fract_prec .* current_parm_est))] on two
+%% consecutive iterations, default = zeros().
+%% Field 'options.max_fract_change': column vector (same length as
+%% 'pin) of maximum fractional step changes in parameter vector.
+%% Fractional change in elements of parameter vector is constrained to
+%% be at most 'options.max_fract_change' between sucessive iterations.
+%% [ie. abs(chg(i))=abs(min([chg(i)
+%% options.max_fract_change(i)*current param estimate])).], default =
+%% Inf*ones().
+%% Field 'options.inequc': cell-array containing up to four entries,
+%% two entries for linear inequality constraints and/or one or two
+%% entries for general inequality constraints. Initial parameters
+%% must satisfy these constraints. Either linear or general
+%% constraints may be the first entries, but the two entries for
+%% linear constraints must be adjacent and, if two entries are given
+%% for general constraints, they also must be adjacent. The two
+%% entries for linear constraints are a matrix (say m) and a vector
+%% (say v), specifying linear inequality constraints of the form
+%% `m.' * parameters + v >= 0'. If the constraints are just bounds,
+%% it is suggested to specify them in 'options.bounds' instead,
+%% since then some sanity tests are performed, and since the
+%% function 'dfdp.m' is guarantied not to violate constraints during
+%% determination of the numeric gradient only for those constraints
+%% specified as 'bounds' (possibly with violations due to a certain
+%% inaccuracy, however, except if no constraints except bounds are
+%% specified). The first entry for general constraints must be a
+%% differentiable vector valued function (say h), specifying general
+%% inequality constraints of the form `h (p[, idx]) >= 0'; p is the
+%% column vector of optimized paraters and the optional argument idx
+%% is a logical index. h has to return the values of all constraints
+%% if idx is not given, and has to return only the indexed
+%% constraints if idx is given (so computation of the other
+%% constraints can be spared). If a second entry for general
+%% constraints is given, it must be a function (say dh) which
+%% returnes a matrix whos rows contain the gradients of the
+%% constraint function h with respect to the optimized parameters.
+%% It has the form jac_h = dh (vh, p, dp, h, idx[, bounds]); p is
+%% the column vector of optimized parameters, and idx is a logical
+%% index --- only the rows indexed by idx must be returned (so
+%% computation of the others can be spared). The other arguments of
+%% dh are for the case that dh computes numerical gradients: vh is
+%% the column vector of the current values of the constraint
+%% function h, with idx already applied. h is a function h (p) to
+%% compute the values of the constraints for parameters p, it will
+%% return only the values indexed by idx. dp is a suggestion for
+%% relative step width, having the same value as the argument 'dp'
+%% of leasqr above. If bounds were specified to leasqr, they are
+%% provided in the argument bounds of dh, to enable their
+%% consideration in determination of numerical gradients. If dh is
+%% not specified to leasqr, numerical gradients are computed in the
+%% same way as with 'dfdp.m' (see above). If some constraints are
+%% linear, they should be specified as linear constraints (or
+%% bounds, if applicable) for reasons of performance, even if
+%% general constraints are also specified.
+%% Field 'options.bounds': two-column-matrix, one row for each
+%% parameter in 'pin'. Each row contains a minimal and maximal value
+%% for each parameter. Default: [-Inf, Inf] in each row. If this
+%% field is used with an existing user-side function for 'dFdp'
+%% (see above) the functions interface might have to be changed.
+%% Field 'options.equc': equality constraints, specified the same
+%% way as inequality constraints (see field 'options.inequc').
+%% Initial parameters must satisfy these constraints.
+%% Note that there is possibly a certain inaccuracy in honoring
+%% constraints, except if only bounds are specified.
+%% _Warning_: If constraints (or bounds) are set, returned guesses
+%% of corp, covp, and Z are generally invalid, even if no constraints
+%% are active for the final parameters. If equality constraints are
+%% specified, corp, covp, and Z are not guessed at all.
+%% Field 'options.cpiv': Function for complementary pivot algorithm
+%% for inequality constraints, default: cpiv_bard. No different
+%% function is supplied.
+%%
+%% OUTPUT VARIABLES
+%% f = column vector of values computed: f = F(x,p).
+%% p = column vector trial or final parameters. i.e, the solution.
+%% cvg = scalar: = 1 if convergence, = 0 otherwise.
+%% iter = scalar number of iterations used.
+%% corp = correlation matrix for parameters.
+%% covp = covariance matrix of the parameters.
+%% covr = diag(covariance matrix of the residuals).
+%% stdresid = standardized residuals.
+%% Z = matrix that defines confidence region (see comments in the source).
+%% r2 = coefficient of multiple determination, intercept form.
+%%
+%% Not suitable for non-real residuals.
+%%
+%% References:
+%% Bard, Nonlinear Parameter Estimation, Academic Press, 1974.
+%% Draper and Smith, Applied Regression Analysis, John Wiley and Sons, 1981.
+
+function [f,p,cvg,iter,corp,covp,covr,stdresid,Z,r2]= ...
+ leasqr(x,y,pin,F,stol,niter,wt,dp,dFdp,options)
+
+ %% The following two blocks of comments are chiefly from the original
+ %% version for Matlab. For later changes the logs of the Octave Forge
+ %% svn repository should also be consulted.
+
+ %% A modified version of Levenberg-Marquardt
+ %% Non-Linear Regression program previously submitted by R.Schrager.
+ %% This version corrects an error in that version and also provides
+ %% an easier to use version with automatic numerical calculation of
+ %% the Jacobian Matrix. In addition, this version calculates statistics
+ %% such as correlation, etc....
+ %%
+ %% Version 3 Notes
+ %% Errors in the original version submitted by Shrager (now called
+ %% version 1) and the improved version of Jutan (now called version 2)
+ %% have been corrected.
+ %% Additional features, statistical tests, and documentation have also been
+ %% included along with an example of usage. BEWARE: Some the the input and
+ %% output arguments were changed from the previous version.
+ %%
+ %% Ray Muzic <rfm2@ds2.uh.cwru.edu>
+ %% Arthur Jutan <jutan@charon.engga.uwo.ca>
+
+ %% Richard I. Shrager (301)-496-1122
+ %% Modified by A.Jutan (519)-679-2111
+ %% Modified by Ray Muzic 14-Jul-1992
+ %% 1) add maxstep feature for limiting changes in parameter estimates
+ %% at each step.
+ %% 2) remove forced columnization of x (x=x(:)) at beginning. x
+ %% could be a matrix with the ith row of containing values of
+ %% the independent variables at the ith observation.
+ %% 3) add verbose option
+ %% 4) add optional return arguments covp, stdresid, chi2
+ %% 5) revise estimates of corp, stdev
+ %% Modified by Ray Muzic 11-Oct-1992
+ %% 1) revise estimate of Vy. remove chi2, add Z as return values
+ %% (later remark: the current code contains no variable Vy)
+ %% Modified by Ray Muzic 7-Jan-1994
+ %% 1) Replace ones(x) with a construct that is compatible with versions
+ %% newer and older than v 4.1.
+ %% 2) Added global declaration of verbose (needed for newer than v4.x)
+ %% 3) Replace return value var, the variance of the residuals
+ %% with covr, the covariance matrix of the residuals.
+ %% 4) Introduce options as 10th input argument. Include
+ %% convergence criteria and maxstep in it.
+ %% 5) Correct calculation of xtx which affects coveraince estimate.
+ %% 6) Eliminate stdev (estimate of standard deviation of
+ %% parameter estimates) from the return values. The covp is a
+ %% much more meaningful expression of precision because it
+ %% specifies a confidence region in contrast to a confidence
+ %% interval.. If needed, however, stdev may be calculated as
+ %% stdev=sqrt(diag(covp)).
+ %% 7) Change the order of the return values to a more logical order.
+ %% 8) Change to more efficent algorithm of Bard for selecting epsL.
+ %% 9) Tighten up memory usage by making use of sparse matrices (if
+ %% MATLAB version >= 4.0) in computation of covp, corp, stdresid.
+ %% Modified by Francesco Potortì
+ %% for use in Octave
+ %% Added linear inequality constraints with quadratic programming to
+ %% this file and special case bounds to this file and to dfdp.m
+ %% (24-Feb-2010) and later also general inequality constraints
+ %% (12-Apr-2010) (Reference: Bard, Y., 'An eclectic approach to
+ %% nonlinear programming', Proc. ANU Sem. Optimization, Canberra,
+ %% Austral. Nat. Univ.). Differences from the reference: adaption to
+ %% svd-based algorithm, linesearch or stepwidth adaptions to ensure
+ %% decrease in objective function omitted to rather start a new
+ %% overall cycle with a new epsL, some performance gains from linear
+ %% constraints even if general constraints are specified. Equality
+ %% constraints also implemented. Olaf Till
+ %% Now split into files leasqr.m and __lm_svd__.m.
+
+ %% needed for some anonymous functions
+ if (exist ('ifelse') ~= 5)
+ ifelse = @ scalar_ifelse;
+ end
+
+ __plot_cmds__ (); % flag persistent variables invalid
+
+ global verbose;
+
+ %% argument processing
+ %%
+
+ if (nargin > 8)
+ if (ischar (dFdp))
+ dfdp = str2func (dFdp);
+ else
+ dfdp = dFdp;
+ end
+ end
+
+ if (nargin <= 7) dp=.001*(pin*0+1); end %DT
+ if (nargin <= 6) wt = ones (size (y)); end % SMB modification
+ if (nargin <= 5) niter = []; end
+ if (nargin == 4) stol=.0001; end
+ if (ischar (F)) F = str2func (F); end
+ %%
+
+ if (any (size (y) ~= size (wt)))
+ error ('dimensions of observations and weights do not match');
+ end
+ wtl = wt(:);
+ pin=pin(:); dp=dp(:); %change all vectors to columns
+ [rows_y, cols_y] = size (y);
+ m = rows_y * cols_y; n=length(pin);
+ f_pin = F (x, pin);
+ if (any (size (f_pin) ~= size (y)))
+ error ('dimensions of returned values of model function and of observations do not match');
+ end
+ f_pin = y - f_pin;
+
+ dFdp = @ (p, dfdp_hook) - dfdp (x, y(:) - dfdp_hook.f, p, dp, F);
+
+ %% processing of 'options'
+ pprec = zeros (n, 1);
+ maxstep = Inf * ones (n, 1);
+ have_gencstr = false; % no general constraints
+ have_genecstr = false; % no general equality constraints
+ n_gencstr = 0;
+ mc = zeros (n, 0);
+ vc = zeros (0, 1); rv = 0;
+ emc = zeros (n, 0);
+ evc = zeros (0, 1); erv = 0;
+ bounds = cat (2, -Inf * ones (n, 1), Inf * ones (n, 1));
+ pin_cstr.inequ.lin_except_bounds = [];;
+ pin_cstr.inequ.gen = [];;
+ pin_cstr.equ.lin = [];;
+ pin_cstr.equ.gen = [];;
+ dfdp_bounds = {};
+ cpiv = @ cpiv_bard;
+ eq_idx = []; % numerical index for equality constraints in all
+ % constraints, later converted to
+ % logical index
+ if (nargin > 9)
+ if (ismatrix (options)) % backwards compatibility
+ tp = options;
+ options = struct ('fract_prec', tp(:, 1));
+ if (columns (tp) > 1)
+ options.max_fract_change = tp(:, 2);
+ end
+ end
+ if (isfield (options, 'cpiv') && ~isempty (options.cpiv))
+ %% As yet there is only one cpiv function distributed with leasqr,
+ %% but this may change; the algorithm of cpiv_bard is said to be
+ %% relatively fast, but may have disadvantages.
+ if (ischar (options.cpiv))
+ cpiv = str2func (options.cpiv);
+ else
+ cpiv = options.cpiv;
+ end
+ end
+ if (isfield (options, 'fract_prec'))
+ pprec = options.fract_prec;
+ if (any (size (pprec) ~= [n, 1]))
+ error ('fractional precisions: wrong dimensions');
+ end
+ end
+ if (isfield (options, 'max_fract_change'))
+ maxstep = options.max_fract_change;
+ if (any (size (maxstep) ~= [n, 1]))
+ error ('maximum fractional step changes: wrong dimensions');
+ end
+ end
+ if (isfield (options, 'inequc'))
+ inequc = options.inequc;
+ if (ismatrix (inequc{1}))
+ mc = inequc{1};
+ vc = inequc{2};
+ if (length (inequc) > 2)
+ have_gencstr = true;
+ f_gencstr = inequc{3};
+ if (length (inequc) > 3)
+ df_gencstr = inequc{4};
+ else
+ df_gencstr = @ dcdp;
+ end
+ end
+ else
+ lid = 0; % no linear constraints
+ have_gencstr = true;
+ f_gencstr = inequc{1};
+ if (length (inequc) > 1)
+ if (ismatrix (inequc{2}))
+ lid = 2;
+ df_gencstr = @ dcdp;
+ else
+ df_gencstr = inequc{2};
+ if (length (inequc) > 2)
+ lid = 3;
+ end
+ end
+ else
+ df_gencstr = @ dcdp;
+ end
+ if (lid)
+ mc = inequc{lid};
+ vc = inequc{lid + 1};
+ end
+ end
+ if (have_gencstr)
+ if (ischar (f_gencstr))
+ f_gencstr = str2func (f_gencstr);
+ end
+ tp = f_gencstr (pin);
+ n_gencstr = length (tp);
+ f_gencstr = @ (p, idx) tf_gencstr (p, idx, f_gencstr);
+ if (ischar (df_gencstr))
+ df_gencstr = str2func (df_gencstr);
+ end
+ if (strcmp (func2str (df_gencstr), 'dcdp'))
+ df_gencstr = @ (f, p, dp, idx, db) ...
+ df_gencstr (f(idx), p, dp, ...
+ @ (tp) f_gencstr (tp, idx), db{:});
+ else
+ df_gencstr = @ (f, p, dp, idx, db) ...
+ df_gencstr (f(idx), p, dp, ...
+ @ (tp) f_gencstr (tp, idx), idx, db{:});
+ end
+ end
+ [rm, cm] = size (mc);
+ [rv, cv] = size (vc);
+ if (rm ~= n || cm ~= rv || cv ~= 1)
+ error ('linear inequality constraints: wrong dimensions');
+ end
+ pin_cstr.inequ.lin_except_bounds = mc.' * pin + vc;
+ if (have_gencstr)
+ pin_cstr.inequ.gen = tp;
+ end
+ end
+ if (isfield (options, 'equc'))
+ equc = options.equc;
+ if (ismatrix (equc{1}))
+ emc = equc{1};
+ evc = equc{2};
+ if (length (equc) > 2)
+ have_genecstr = true;
+ f_genecstr = equc{3};
+ if (length (equc) > 3)
+ df_genecstr = equc{4};
+ else
+ df_genecstr = @ dcdp;
+ end
+ end
+ else
+ lid = 0; % no linear constraints
+ have_genecstr = true;
+ f_genecstr = equc{1};
+ if (length (equc) > 1)
+ if (ismatrix (equc{2}))
+ lid = 2;
+ df_genecstr = @ dcdp;
+ else
+ df_genecstr = equc{2};
+ if (length (equc) > 2)
+ lid = 3;
+ end
+ end
+ else
+ df_genecstr = @ dcdp;
+ end
+ if (lid)
+ emc = equc{lid};
+ evc = equc{lid + 1};
+ end
+ end
+ if (have_genecstr)
+ if (ischar (f_genecstr))
+ f_genecstr = str2func (f_genecstr);
+ end
+ tp = f_genecstr (pin);
+ n_genecstr = length (tp);
+ f_genecstr = @ (p, idx) tf_gencstr (p, idx, f_genecstr);
+ if (ischar (df_genecstr))
+ df_genecstr = str2func (df_genecstr);
+ end
+ if (strcmp (func2str (df_genecstr), 'dcdp'))
+ df_genecstr = @ (f, p, dp, idx, db) ...
+ df_genecstr (f, p, dp, ...
+ @ (tp) f_genecstr (tp, idx), db{:});
+ else
+ df_genecstr = @ (f, p, dp, idx, db) ...
+ df_genecstr (f, p, dp, ...
+ @ (tp) f_genecstr (tp, idx), idx, db{:});
+ end
+ end
+ [erm, ecm] = size (emc);
+ [erv, ecv] = size (evc);
+ if (erm ~= n || ecm ~= erv || ecv ~= 1)
+ error ('linear equality constraints: wrong dimensions');
+ end
+ pin_cstr.equ.lin = emc.' * pin + evc;
+ if (have_genecstr)
+ pin_cstr.equ.gen = tp;
+ end
+ end
+ if (isfield (options, 'bounds'))
+ bounds = options.bounds;
+ if (any (size (bounds) ~= [n, 2]))
+ error ('bounds: wrong dimensions');
+ end
+ idx = bounds(:, 1) > bounds(:, 2);
+ tp = bounds(idx, 2);
+ bounds(idx, 2) = bounds(idx, 1);
+ bounds(idx, 1) = tp;
+ %% It is possible to take this decision here, since this frontend
+ %% is used only with one certain backend. The backend will check
+ %% this again; but it will not reference 'dp' in its message,
+ %% thats why the additional check here.
+ idx = bounds(:, 1) == bounds(:, 2);
+ if (any (idx))
+ warning ('leasqr:constraints', 'lower and upper bounds identical for some parameters, setting the respective elements of dp to zero');
+ dp(idx) = 0;
+ end
+ %%
+ tp = eye (n);
+ lidx = ~isinf (bounds(:, 1));
+ uidx = ~isinf (bounds(:, 2));
+ mc = cat (2, mc, tp(:, lidx), - tp(:, uidx));
+ vc = cat (1, vc, - bounds(lidx, 1), bounds(uidx, 2));
+ [rm, cm] = size (mc);
+ [rv, cv] = size (vc);
+ dfdp_bounds = {bounds};
+ dFdp = @ (p, dfdp_hook) - dfdp (x, y(:) - dfdp_hook.f, p, dp, ...
+ F, bounds);
+ end
+ %% concatenate inequality and equality constraint functions, mc, and
+ %% vc; update eq_idx, rv, n_gencstr, have_gencstr
+ if (erv > 0)
+ mc = cat (2, mc, emc);
+ vc = cat (1, vc, evc);
+ eq_idx = rv + 1 : rv + erv;
+ rv = rv + erv;
+ end
+ if (have_genecstr)
+ eq_idx = cat (2, eq_idx, ...
+ rv + n_gencstr + 1 : rv + n_gencstr + n_genecstr);
+ nidxi = 1 : n_gencstr;
+ nidxe = n_gencstr + 1 : n_gencstr + n_genecstr;
+ n_gencstr = n_gencstr + n_genecstr;
+ if (have_gencstr)
+ f_gencstr = @ (p, idx) cat (1, ...
+ f_gencstr (p, idx(nidxi)), ...
+ f_genecstr (p, idx(nidxe)));
+ df_gencstr = @ (f, p, dp, idx, db) ...
+ cat (1, ...
+ df_gencstr (f(nidxi), p, dp, idx(nidxi), db), ...
+ df_genecstr (f(nidxe), p, dp, idx(nidxe), db));
+ else
+ f_gencstr = f_genecstr;
+ df_gencstr = df_genecstr;
+ have_gencstr = true;
+ end
+ end
+ end
+ if (have_gencstr)
+ nidxl = 1:rv;
+ nidxh = rv+1:rv+n_gencstr;
+ f_cstr = @ (p, idx) ...
+ cat (1, mc(:, idx(nidxl)).' * p + vc(idx(nidxl), 1), ...
+ f_gencstr (p, idx(nidxh)));
+ %% in the case of this interface, diffp is already zero at fixed;
+ %% also in this special case, dfdp_bounds can be filled in directly
+ %% --- otherwise it would be a field of hook in the called function
+ df_cstr = @ (p, idx, dfdp_hook) ...
+ cat (1, mc(:, idx(nidxl)).', ...
+ df_gencstr (dfdp_hook.f(nidxh), p, dp, ...
+ idx(nidxh), ...
+ dfdp_bounds));
+ else
+ f_cstr = @ (p, idx) mc(:, idx).' * p + vc(idx, 1);
+ df_cstr = @ (p, idx, dfdp_hook) mc(:, idx).';
+ end
+
+
+
+ %% in a general interface, check for all(fixed) here
+
+ %% passed constraints
+ hook.mc = mc; % matrix of linear constraints
+ hook.vc = vc; % vector of linear constraints
+ hook.f_cstr = f_cstr; % function of all constraints
+ hook.df_cstr = df_cstr; % function of derivatives of all constraints
+ hook.n_gencstr = n_gencstr; % number of non-linear constraints
+ hook.eq_idx = false (size (vc, 1) + n_gencstr, 1);
+ hook.eq_idx(eq_idx) = true; % logical index of equality constraints in
+ % all constraints
+ hook.lbound = bounds(:, 1); % bounds, subset of linear inequality
+ % constraints in mc and vc
+ hook.ubound = bounds(:, 2);
+
+ %% passed values of constraints for initial parameters
+ hook.pin_cstr = pin_cstr;
+
+ %% passed derivative of model function
+ hook.dfdp = dFdp;
+
+ %% passed function for complementary pivoting
+ hook.cpiv = cpiv;
+
+ %% passed value of residual function for initial parameters
+ hook.f_pin = f_pin;
+
+ %% passed options
+ hook.max_fract_change = maxstep;
+ hook.fract_prec = pprec;
+ hook.TolFun = stol;
+ hook.MaxIter = niter;
+ hook.weights = wt;
+ hook.fixed = dp == 0;
+ if (verbose)
+ hook.Display = 'iter';
+ __plot_cmds__ = @ __plot_cmds__; # for bug #31484 (Octave <= 3.2.4)
+ hook.plot_cmd = @ (f) __plot_cmds__ (x, y, y - f);
+ else
+ hook.Display = 'off';
+ end
+
+ %% only preliminary, for testing
+ hook.testing = false;
+ hook.new_s = false;
+ if (exist ('options'))
+ if (isfield (options, 'testing'))
+ hook.testing = options.testing;
+ end
+ if (isfield (options, 'new_s'))
+ hook.new_s = options.new_s;
+ end
+ end
+
+ [p, resid, cvg, outp] = __lm_svd__ (@ (p) y - F (x, p), pin, hook);
+ f = y - resid;
+ iter = outp.niter;
+ cvg = cvg > 0;
+
+ if (~cvg) disp(' CONVERGENCE NOT ACHIEVED! '); end
+
+ if (~(verbose || nargout > 4)) return; end
+
+ yl = y(:);
+ f = f(:);
+ %% CALC VARIANCE COV MATRIX AND CORRELATION MATRIX OF PARAMETERS
+ %% re-evaluate the Jacobian at optimal values
+ jac = dFdp (p, struct ('f', f));
+ msk = ~hook.fixed;
+ n = sum(msk); % reduce n to equal number of estimated parameters
+ jac = jac(:, msk); % use only fitted parameters
+
+ %% following section is Ray Muzic's estimate for covariance and correlation
+ %% assuming covariance of data is a diagonal matrix proportional to
+ %% diag(1/wt.^2).
+ %% cov matrix of data est. from Bard Eq. 7-5-13, and Row 1 Table 5.1
+
+ tp = wtl.^2;
+ if (exist('sparse')) % save memory
+ Q = sparse (1:m, 1:m, 1 ./ tp);
+ Qinv = sparse (1:m, 1:m, tp);
+ else
+ Q = diag (ones (m, 1) ./ tp);
+ Qinv = diag (tp);
+ end
+ resid = resid(:); % un-weighted residuals
+ if (~isreal (resid)) error ('residuals are not real'); end
+ tp = resid.' * Qinv * resid;
+ covr = (tp / m) * Q; %covariance of residuals
+
+ %% Matlab compatibility and avoiding recomputation make the following
+ %% logic clumsy.
+ compute = 1;
+ if (m <= n || any (eq_idx))
+ compute = 0;
+ else
+ Qinv = ((m - n) / tp) * Qinv;
+ %% simplified Eq. 7-5-13, Bard; cov of parm est, inverse; outer
+ %% parantheses contain inverse of guessed covariance matrix of data
+ covpinv = jac.' * Qinv * jac;
+ if (exist ('rcond') && rcond (covpinv) <= eps)
+ compute = 0;
+ elseif (rank (covpinv) < n)
+ %% above test is not equivalent to 'rcond' and may unnecessarily
+ %% reject some matrices
+ compute = 0;
+ end
+ end
+ if (compute)
+ covp = inv (covpinv);
+ d=sqrt(diag(covp));
+ corp = covp ./ (d * d.');
+ else
+ covp = NA * ones (n);
+ corp = covp;
+ end
+
+ if (exist('sparse'))
+ covr=spdiags(covr,0);
+ else
+ covr=diag(covr); % convert returned values to
+ % compact storage
+ end
+ covr = reshape (covr, rows_y, cols_y);
+ stdresid = resid .* abs (wtl) / sqrt (tp / m); % equivalent to resid ./
+ % sqrt (covr)
+ stdresid = reshape (stdresid, rows_y, cols_y);
+
+ if (~(verbose || nargout > 8)) return; end
+
+ if (m > n && ~any (eq_idx))
+ Z = ((m - n) / (n * resid.' * Qinv * resid)) * covpinv;
+ else
+ Z = NA * ones (n);
+ end
+
+%%% alt. est. of cov. mat. of parm.:(Delforge, Circulation, 82:1494-1504, 1990
+ %%disp('Alternate estimate of cov. of param. est.')
+ %%acovp=resid'*Qinv*resid/(m-n)*inv(jac'*Qinv*jac);
+
+ if (~(verbose || nargout > 9)) return; end
+
+ %%Calculate R^2, intercept form
+ %%
+ tp = sumsq (yl - mean (yl));
+ if (tp > 0)
+ r2 = 1 - sumsq (resid) / tp;
+ else
+ r2 = NA;
+ end
+
+ %% if someone has asked for it, let them have it
+ %%
+ if (verbose)
+ __plot_cmds__ (x, y, f);
+ disp(' Least Squares Estimates of Parameters')
+ disp(p.')
+ disp(' Correlation matrix of parameters estimated')
+ disp(corp)
+ disp(' Covariance matrix of Residuals' )
+ disp(covr)
+ disp(' Correlation Coefficient R^2')
+ disp(r2)
+ fprintf(' 95%% conf region: F(0.05)(%.0f,%.0f)>= delta_pvec.%s*Z*delta_pvec\n', n, m - n, char (39)); % works with ' and '
+ Z
+ %% runs test according to Bard. p 201.
+ n1 = sum (resid > 0);
+ n2 = sum (resid < 0);
+ nrun=sum(abs(diff(resid > 0)))+1;
+ if ((n1 > 10) && (n2 > 10)) % sufficent data for test?
+ zed=(nrun-(2*n1*n2/(n1+n2)+1)+0.5)/(2*n1*n2*(2*n1*n2-n1-n2)...
+ /((n1+n2)^2*(n1+n2-1)));
+ if (zed < 0)
+ prob = erfc(-zed/sqrt(2))/2*100;
+ disp([num2str(prob),'% chance of fewer than ',num2str(nrun),' runs.']);
+ else
+ prob = erfc(zed/sqrt(2))/2*100;
+ disp([num2str(prob),'% chance of greater than ',num2str(nrun),' runs.']);
+ end
+ end
+ end
+
+function ret = tf_gencstr (p, idx, f)
+
+ %% necessary since user function f_gencstr might return [] or a row
+ %% vector
+
+ ret = f (p, idx);
+ if (isempty (ret))
+ ret = zeros (0, 1);
+ elseif (size (ret, 2) > 1)
+ ret = ret(:);
+ end
+
+function fval = scalar_ifelse (cond, tval, fval)
+
+ %% needed for some anonymous functions, builtin ifelse only available
+ %% in Octave > 3.2; we need only the scalar case here
+
+ if (cond)
+ fval = tval;
+ end
+
+%!demo
+%! % Define functions
+%! leasqrfunc = @(x, p) p(1) * exp (-p(2) * x);
+%! leasqrdfdp = @(x, f, p, dp, func) [exp(-p(2)*x), -p(1)*x.*exp(-p(2)*x)];
+%!
+%! % generate test data
+%! t = [1:10:100]';
+%! p = [1; 0.1];
+%! data = leasqrfunc (t, p);
+%!
+%! rnd = [0.352509; -0.040607; -1.867061; -1.561283; 1.473191; ...
+%! 0.580767; 0.841805; 1.632203; -0.179254; 0.345208];
+%!
+%! % add noise
+%! % wt1 = 1 /sqrt of variances of data
+%! % 1 / wt1 = sqrt of var = standard deviation
+%! wt1 = (1 + 0 * t) ./ sqrt (data);
+%! data = data + 0.05 * rnd ./ wt1;
+%!
+%! % Note by Thomas Walter <walter@pctc.chemie.uni-erlangen.de>:
+%! %
+%! % Using a step size of 1 to calculate the derivative is WRONG !!!!
+%! % See numerical mathbooks why.
+%! % A derivative calculated from central differences need: s
+%! % step = 0.001...1.0e-8
+%! % And onesided derivative needs:
+%! % step = 1.0e-5...1.0e-8 and may be still wrong
+%!
+%! F = leasqrfunc;
+%! dFdp = leasqrdfdp; % exact derivative
+%! % dFdp = dfdp; % estimated derivative
+%! dp = [0.001; 0.001];
+%! pin = [.8; .05];
+%! stol=0.001; niter=50;
+%! minstep = [0.01; 0.01];
+%! maxstep = [0.8; 0.8];
+%! options = [minstep, maxstep];
+%!
+%! global verbose;
+%! verbose = 1;
+%! [f1, p1, kvg1, iter1, corp1, covp1, covr1, stdresid1, Z1, r21] = ...
+%! leasqr (t, data, pin, F, stol, niter, wt1, dp, dFdp, options);
+
+%!demo
+%! %% Example for linear inequality constraints.
+%! %% model function:
+%! F = @ (x, p) p(1) * exp (p(2) * x);
+%! %% independents and dependents:
+%! x = 1:5;
+%! y = [1, 2, 4, 7, 14];
+%! %% initial values:
+%! init = [.25; .25];
+%! %% other configuration (default values):
+%! tolerance = .0001;
+%! max_iterations = 20;
+%! weights = ones (1, 5);
+%! dp = [.001; .001]; % bidirectional numeric gradient stepsize
+%! dFdp = 'dfdp'; % function for gradient (numerical)
+%!
+%! %% linear constraints, A.' * parametervector + B >= 0
+%! A = [1; -1]; B = 0; % p(1) >= p(2);
+%! options.inequc = {A, B};
+%!
+%! %% start leasqr, be sure that 'verbose' is not set
+%! global verbose; verbose = false;
+%! [f, p, cvg, iter] = ...
+%! leasqr (x, y, init, F, tolerance, max_iterations, ...
+%! weights, dp, dFdp, options)
git://git.creatis.insa-lyon.fr / CreaPhase.git / blobdiff