1 %% Copyright (C) 1992-1994 Richard Shrager
2 %% Copyright (C) 1992-1994 Arthur Jutan <jutan@charon.engga.uwo.ca>
3 %% Copyright (C) 1992-1994 Ray Muzic <rfm2@ds2.uh.cwru.edu>
4 %% Copyright (C) 2010, 2011 Olaf Till <i7tiol@t-online.de>
6 %% This program is free software; you can redistribute it and/or modify it under
7 %% the terms of the GNU General Public License as published by the Free Software
8 %% Foundation; either version 3 of the License, or (at your option) any later
11 %% This program is distributed in the hope that it will be useful, but WITHOUT
12 %% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 %% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16 %% You should have received a copy of the GNU General Public License along with
17 %% this program; if not, see <http://www.gnu.org/licenses/>.
19 %%function [f,p,cvg,iter,corp,covp,covr,stdresid,Z,r2]=
20 %% leasqr(x,y,pin,F,{stol,niter,wt,dp,dFdp,options})
22 %% Levenberg-Marquardt nonlinear regression of f(x,p) to y(x).
25 %% Optional parameters are in braces {}.
26 %% x = vector or matrix of independent variables.
27 %% y = vector or matrix of observed values.
28 %% wt = statistical weights (same dimensions as y). These should be
29 %% set to be proportional to (sqrt of var(y))^-1; (That is, the
30 %% covariance matrix of the data is assumed to be proportional to
31 %% diagonal with diagonal equal to (wt.^2)^-1. The constant of
32 %% proportionality will be estimated.); default = ones( size (y)).
33 %% pin = vec of initial parameters to be adjusted by leasqr.
34 %% dp = fractional increment of p for numerical partial derivatives;
35 %% default = .001*ones(size(pin))
36 %% dp(j) > 0 means central differences on j-th parameter p(j).
37 %% dp(j) < 0 means one-sided differences on j-th parameter p(j).
38 %% dp(j) = 0 holds p(j) fixed i.e. leasqr wont change initial guess: pin(j)
39 %% F = name of function in quotes or function handle; the function
40 %% shall be of the form y=f(x,p), with y, x, p of the form y, x, pin
41 %% as described above.
42 %% dFdp = name of partial derivative function in quotes or function
43 %% handle; default is 'dfdp', a slow but general partial derivatives
44 %% function; the function shall be of the form
45 %% prt=dfdp(x,f,p,dp,F[,bounds]). For backwards compatibility, the
46 %% function will only be called with an extra 'bounds' argument if the
47 %% 'bounds' option is explicitely specified to leasqr (see dfdp.m).
48 %% stol = scalar tolerance on fractional improvement in scalar sum of
49 %% squares = sum((wt.*(y-f))^2); default stol = .0001;
50 %% niter = scalar maximum number of iterations; default = 20;
51 %% options = structure, currently recognized fields are 'fract_prec',
52 %% 'max_fract_change', 'inequc', 'bounds', and 'equc'. For backwards
53 %% compatibility, 'options' can also be a matrix whose first and
54 %% second column contains the values of 'fract_prec' and
55 %% 'max_fract_change', respectively.
56 %% Field 'options.fract_prec': column vector (same length as 'pin')
57 %% of desired fractional precisions in parameter estimates.
58 %% Iterations are terminated if change in parameter vector (chg)
59 %% relative to current parameter estimate is less than their
60 %% corresponding elements in 'options.fract_prec' [ie. all (abs
61 %% (chg) < abs (options.fract_prec .* current_parm_est))] on two
62 %% consecutive iterations, default = zeros().
63 %% Field 'options.max_fract_change': column vector (same length as
64 %% 'pin) of maximum fractional step changes in parameter vector.
65 %% Fractional change in elements of parameter vector is constrained to
66 %% be at most 'options.max_fract_change' between sucessive iterations.
67 %% [ie. abs(chg(i))=abs(min([chg(i)
68 %% options.max_fract_change(i)*current param estimate])).], default =
70 %% Field 'options.inequc': cell-array containing up to four entries,
71 %% two entries for linear inequality constraints and/or one or two
72 %% entries for general inequality constraints. Initial parameters
73 %% must satisfy these constraints. Either linear or general
74 %% constraints may be the first entries, but the two entries for
75 %% linear constraints must be adjacent and, if two entries are given
76 %% for general constraints, they also must be adjacent. The two
77 %% entries for linear constraints are a matrix (say m) and a vector
78 %% (say v), specifying linear inequality constraints of the form
79 %% `m.' * parameters + v >= 0'. If the constraints are just bounds,
80 %% it is suggested to specify them in 'options.bounds' instead,
81 %% since then some sanity tests are performed, and since the
82 %% function 'dfdp.m' is guarantied not to violate constraints during
83 %% determination of the numeric gradient only for those constraints
84 %% specified as 'bounds' (possibly with violations due to a certain
85 %% inaccuracy, however, except if no constraints except bounds are
86 %% specified). The first entry for general constraints must be a
87 %% differentiable vector valued function (say h), specifying general
88 %% inequality constraints of the form `h (p[, idx]) >= 0'; p is the
89 %% column vector of optimized paraters and the optional argument idx
90 %% is a logical index. h has to return the values of all constraints
91 %% if idx is not given, and has to return only the indexed
92 %% constraints if idx is given (so computation of the other
93 %% constraints can be spared). If a second entry for general
94 %% constraints is given, it must be a function (say dh) which
95 %% returnes a matrix whos rows contain the gradients of the
96 %% constraint function h with respect to the optimized parameters.
97 %% It has the form jac_h = dh (vh, p, dp, h, idx[, bounds]); p is
98 %% the column vector of optimized parameters, and idx is a logical
99 %% index --- only the rows indexed by idx must be returned (so
100 %% computation of the others can be spared). The other arguments of
101 %% dh are for the case that dh computes numerical gradients: vh is
102 %% the column vector of the current values of the constraint
103 %% function h, with idx already applied. h is a function h (p) to
104 %% compute the values of the constraints for parameters p, it will
105 %% return only the values indexed by idx. dp is a suggestion for
106 %% relative step width, having the same value as the argument 'dp'
107 %% of leasqr above. If bounds were specified to leasqr, they are
108 %% provided in the argument bounds of dh, to enable their
109 %% consideration in determination of numerical gradients. If dh is
110 %% not specified to leasqr, numerical gradients are computed in the
111 %% same way as with 'dfdp.m' (see above). If some constraints are
112 %% linear, they should be specified as linear constraints (or
113 %% bounds, if applicable) for reasons of performance, even if
114 %% general constraints are also specified.
115 %% Field 'options.bounds': two-column-matrix, one row for each
116 %% parameter in 'pin'. Each row contains a minimal and maximal value
117 %% for each parameter. Default: [-Inf, Inf] in each row. If this
118 %% field is used with an existing user-side function for 'dFdp'
119 %% (see above) the functions interface might have to be changed.
120 %% Field 'options.equc': equality constraints, specified the same
121 %% way as inequality constraints (see field 'options.inequc').
122 %% Initial parameters must satisfy these constraints.
123 %% Note that there is possibly a certain inaccuracy in honoring
124 %% constraints, except if only bounds are specified.
125 %% _Warning_: If constraints (or bounds) are set, returned guesses
126 %% of corp, covp, and Z are generally invalid, even if no constraints
127 %% are active for the final parameters. If equality constraints are
128 %% specified, corp, covp, and Z are not guessed at all.
129 %% Field 'options.cpiv': Function for complementary pivot algorithm
130 %% for inequality constraints, default: cpiv_bard. No different
131 %% function is supplied.
134 %% f = column vector of values computed: f = F(x,p).
135 %% p = column vector trial or final parameters. i.e, the solution.
136 %% cvg = scalar: = 1 if convergence, = 0 otherwise.
137 %% iter = scalar number of iterations used.
138 %% corp = correlation matrix for parameters.
139 %% covp = covariance matrix of the parameters.
140 %% covr = diag(covariance matrix of the residuals).
141 %% stdresid = standardized residuals.
142 %% Z = matrix that defines confidence region (see comments in the source).
143 %% r2 = coefficient of multiple determination, intercept form.
145 %% Not suitable for non-real residuals.
148 %% Bard, Nonlinear Parameter Estimation, Academic Press, 1974.
149 %% Draper and Smith, Applied Regression Analysis, John Wiley and Sons, 1981.
151 function [f,p,cvg,iter,corp,covp,covr,stdresid,Z,r2]= ...
152 leasqr(x,y,pin,F,stol,niter,wt,dp,dFdp,options)
154 %% The following two blocks of comments are chiefly from the original
155 %% version for Matlab. For later changes the logs of the Octave Forge
156 %% svn repository should also be consulted.
158 %% A modified version of Levenberg-Marquardt
159 %% Non-Linear Regression program previously submitted by R.Schrager.
160 %% This version corrects an error in that version and also provides
161 %% an easier to use version with automatic numerical calculation of
162 %% the Jacobian Matrix. In addition, this version calculates statistics
163 %% such as correlation, etc....
166 %% Errors in the original version submitted by Shrager (now called
167 %% version 1) and the improved version of Jutan (now called version 2)
168 %% have been corrected.
169 %% Additional features, statistical tests, and documentation have also been
170 %% included along with an example of usage. BEWARE: Some the the input and
171 %% output arguments were changed from the previous version.
173 %% Ray Muzic <rfm2@ds2.uh.cwru.edu>
174 %% Arthur Jutan <jutan@charon.engga.uwo.ca>
176 %% Richard I. Shrager (301)-496-1122
177 %% Modified by A.Jutan (519)-679-2111
178 %% Modified by Ray Muzic 14-Jul-1992
179 %% 1) add maxstep feature for limiting changes in parameter estimates
181 %% 2) remove forced columnization of x (x=x(:)) at beginning. x
182 %% could be a matrix with the ith row of containing values of
183 %% the independent variables at the ith observation.
184 %% 3) add verbose option
185 %% 4) add optional return arguments covp, stdresid, chi2
186 %% 5) revise estimates of corp, stdev
187 %% Modified by Ray Muzic 11-Oct-1992
188 %% 1) revise estimate of Vy. remove chi2, add Z as return values
189 %% (later remark: the current code contains no variable Vy)
190 %% Modified by Ray Muzic 7-Jan-1994
191 %% 1) Replace ones(x) with a construct that is compatible with versions
192 %% newer and older than v 4.1.
193 %% 2) Added global declaration of verbose (needed for newer than v4.x)
194 %% 3) Replace return value var, the variance of the residuals
195 %% with covr, the covariance matrix of the residuals.
196 %% 4) Introduce options as 10th input argument. Include
197 %% convergence criteria and maxstep in it.
198 %% 5) Correct calculation of xtx which affects coveraince estimate.
199 %% 6) Eliminate stdev (estimate of standard deviation of
200 %% parameter estimates) from the return values. The covp is a
201 %% much more meaningful expression of precision because it
202 %% specifies a confidence region in contrast to a confidence
203 %% interval.. If needed, however, stdev may be calculated as
204 %% stdev=sqrt(diag(covp)).
205 %% 7) Change the order of the return values to a more logical order.
206 %% 8) Change to more efficent algorithm of Bard for selecting epsL.
207 %% 9) Tighten up memory usage by making use of sparse matrices (if
208 %% MATLAB version >= 4.0) in computation of covp, corp, stdresid.
209 %% Modified by Francesco Potortì
211 %% Added linear inequality constraints with quadratic programming to
212 %% this file and special case bounds to this file and to dfdp.m
213 %% (24-Feb-2010) and later also general inequality constraints
214 %% (12-Apr-2010) (Reference: Bard, Y., 'An eclectic approach to
215 %% nonlinear programming', Proc. ANU Sem. Optimization, Canberra,
216 %% Austral. Nat. Univ.). Differences from the reference: adaption to
217 %% svd-based algorithm, linesearch or stepwidth adaptions to ensure
218 %% decrease in objective function omitted to rather start a new
219 %% overall cycle with a new epsL, some performance gains from linear
220 %% constraints even if general constraints are specified. Equality
221 %% constraints also implemented. Olaf Till
222 %% Now split into files leasqr.m and __lm_svd__.m.
224 %% needed for some anonymous functions
225 if (exist ('ifelse') ~= 5)
226 ifelse = @ scalar_ifelse;
229 __plot_cmds__ (); % flag persistent variables invalid
233 %% argument processing
238 dfdp = str2func (dFdp);
244 if (nargin <= 7) dp=.001*(pin*0+1); end %DT
245 if (nargin <= 6) wt = ones (size (y)); end % SMB modification
246 if (nargin <= 5) niter = []; end
247 if (nargin == 4) stol=.0001; end
248 if (ischar (F)) F = str2func (F); end
251 if (any (size (y) ~= size (wt)))
252 error ('dimensions of observations and weights do not match');
255 pin=pin(:); dp=dp(:); %change all vectors to columns
256 [rows_y, cols_y] = size (y);
257 m = rows_y * cols_y; n=length(pin);
259 if (any (size (f_pin) ~= size (y)))
260 error ('dimensions of returned values of model function and of observations do not match');
264 dFdp = @ (p, dfdp_hook) - dfdp (x, y(:) - dfdp_hook.f, p, dp, F);
266 %% processing of 'options'
267 pprec = zeros (n, 1);
268 maxstep = Inf * ones (n, 1);
269 have_gencstr = false; % no general constraints
270 have_genecstr = false; % no general equality constraints
273 vc = zeros (0, 1); rv = 0;
275 evc = zeros (0, 1); erv = 0;
276 bounds = cat (2, -Inf * ones (n, 1), Inf * ones (n, 1));
277 pin_cstr.inequ.lin_except_bounds = [];;
278 pin_cstr.inequ.gen = [];;
279 pin_cstr.equ.lin = [];;
280 pin_cstr.equ.gen = [];;
283 eq_idx = []; % numerical index for equality constraints in all
284 % constraints, later converted to
287 if (ismatrix (options)) % backwards compatibility
289 options = struct ('fract_prec', tp(:, 1));
290 if (columns (tp) > 1)
291 options.max_fract_change = tp(:, 2);
294 if (isfield (options, 'cpiv') && ~isempty (options.cpiv))
295 %% As yet there is only one cpiv function distributed with leasqr,
296 %% but this may change; the algorithm of cpiv_bard is said to be
297 %% relatively fast, but may have disadvantages.
298 if (ischar (options.cpiv))
299 cpiv = str2func (options.cpiv);
304 if (isfield (options, 'fract_prec'))
305 pprec = options.fract_prec;
306 if (any (size (pprec) ~= [n, 1]))
307 error ('fractional precisions: wrong dimensions');
310 if (isfield (options, 'max_fract_change'))
311 maxstep = options.max_fract_change;
312 if (any (size (maxstep) ~= [n, 1]))
313 error ('maximum fractional step changes: wrong dimensions');
316 if (isfield (options, 'inequc'))
317 inequc = options.inequc;
318 if (ismatrix (inequc{1}))
321 if (length (inequc) > 2)
323 f_gencstr = inequc{3};
324 if (length (inequc) > 3)
325 df_gencstr = inequc{4};
331 lid = 0; % no linear constraints
333 f_gencstr = inequc{1};
334 if (length (inequc) > 1)
335 if (ismatrix (inequc{2}))
339 df_gencstr = inequc{2};
340 if (length (inequc) > 2)
349 vc = inequc{lid + 1};
353 if (ischar (f_gencstr))
354 f_gencstr = str2func (f_gencstr);
356 tp = f_gencstr (pin);
357 n_gencstr = length (tp);
358 f_gencstr = @ (p, idx) tf_gencstr (p, idx, f_gencstr);
359 if (ischar (df_gencstr))
360 df_gencstr = str2func (df_gencstr);
362 if (strcmp (func2str (df_gencstr), 'dcdp'))
363 df_gencstr = @ (f, p, dp, idx, db) ...
364 df_gencstr (f(idx), p, dp, ...
365 @ (tp) f_gencstr (tp, idx), db{:});
367 df_gencstr = @ (f, p, dp, idx, db) ...
368 df_gencstr (f(idx), p, dp, ...
369 @ (tp) f_gencstr (tp, idx), idx, db{:});
372 [rm, cm] = size (mc);
373 [rv, cv] = size (vc);
374 if (rm ~= n || cm ~= rv || cv ~= 1)
375 error ('linear inequality constraints: wrong dimensions');
377 pin_cstr.inequ.lin_except_bounds = mc.' * pin + vc;
379 pin_cstr.inequ.gen = tp;
382 if (isfield (options, 'equc'))
384 if (ismatrix (equc{1}))
387 if (length (equc) > 2)
388 have_genecstr = true;
389 f_genecstr = equc{3};
390 if (length (equc) > 3)
391 df_genecstr = equc{4};
393 df_genecstr = @ dcdp;
397 lid = 0; % no linear constraints
398 have_genecstr = true;
399 f_genecstr = equc{1};
400 if (length (equc) > 1)
401 if (ismatrix (equc{2}))
403 df_genecstr = @ dcdp;
405 df_genecstr = equc{2};
406 if (length (equc) > 2)
411 df_genecstr = @ dcdp;
419 if (ischar (f_genecstr))
420 f_genecstr = str2func (f_genecstr);
422 tp = f_genecstr (pin);
423 n_genecstr = length (tp);
424 f_genecstr = @ (p, idx) tf_gencstr (p, idx, f_genecstr);
425 if (ischar (df_genecstr))
426 df_genecstr = str2func (df_genecstr);
428 if (strcmp (func2str (df_genecstr), 'dcdp'))
429 df_genecstr = @ (f, p, dp, idx, db) ...
430 df_genecstr (f, p, dp, ...
431 @ (tp) f_genecstr (tp, idx), db{:});
433 df_genecstr = @ (f, p, dp, idx, db) ...
434 df_genecstr (f, p, dp, ...
435 @ (tp) f_genecstr (tp, idx), idx, db{:});
438 [erm, ecm] = size (emc);
439 [erv, ecv] = size (evc);
440 if (erm ~= n || ecm ~= erv || ecv ~= 1)
441 error ('linear equality constraints: wrong dimensions');
443 pin_cstr.equ.lin = emc.' * pin + evc;
445 pin_cstr.equ.gen = tp;
448 if (isfield (options, 'bounds'))
449 bounds = options.bounds;
450 if (any (size (bounds) ~= [n, 2]))
451 error ('bounds: wrong dimensions');
453 idx = bounds(:, 1) > bounds(:, 2);
455 bounds(idx, 2) = bounds(idx, 1);
457 %% It is possible to take this decision here, since this frontend
458 %% is used only with one certain backend. The backend will check
459 %% this again; but it will not reference 'dp' in its message,
460 %% thats why the additional check here.
461 idx = bounds(:, 1) == bounds(:, 2);
463 warning ('leasqr:constraints', 'lower and upper bounds identical for some parameters, setting the respective elements of dp to zero');
468 lidx = ~isinf (bounds(:, 1));
469 uidx = ~isinf (bounds(:, 2));
470 mc = cat (2, mc, tp(:, lidx), - tp(:, uidx));
471 vc = cat (1, vc, - bounds(lidx, 1), bounds(uidx, 2));
472 [rm, cm] = size (mc);
473 [rv, cv] = size (vc);
474 dfdp_bounds = {bounds};
475 dFdp = @ (p, dfdp_hook) - dfdp (x, y(:) - dfdp_hook.f, p, dp, ...
478 %% concatenate inequality and equality constraint functions, mc, and
479 %% vc; update eq_idx, rv, n_gencstr, have_gencstr
481 mc = cat (2, mc, emc);
482 vc = cat (1, vc, evc);
483 eq_idx = rv + 1 : rv + erv;
487 eq_idx = cat (2, eq_idx, ...
488 rv + n_gencstr + 1 : rv + n_gencstr + n_genecstr);
489 nidxi = 1 : n_gencstr;
490 nidxe = n_gencstr + 1 : n_gencstr + n_genecstr;
491 n_gencstr = n_gencstr + n_genecstr;
493 f_gencstr = @ (p, idx) cat (1, ...
494 f_gencstr (p, idx(nidxi)), ...
495 f_genecstr (p, idx(nidxe)));
496 df_gencstr = @ (f, p, dp, idx, db) ...
498 df_gencstr (f(nidxi), p, dp, idx(nidxi), db), ...
499 df_genecstr (f(nidxe), p, dp, idx(nidxe), db));
501 f_gencstr = f_genecstr;
502 df_gencstr = df_genecstr;
509 nidxh = rv+1:rv+n_gencstr;
510 f_cstr = @ (p, idx) ...
511 cat (1, mc(:, idx(nidxl)).' * p + vc(idx(nidxl), 1), ...
512 f_gencstr (p, idx(nidxh)));
513 %% in the case of this interface, diffp is already zero at fixed;
514 %% also in this special case, dfdp_bounds can be filled in directly
515 %% --- otherwise it would be a field of hook in the called function
516 df_cstr = @ (p, idx, dfdp_hook) ...
517 cat (1, mc(:, idx(nidxl)).', ...
518 df_gencstr (dfdp_hook.f(nidxh), p, dp, ...
522 f_cstr = @ (p, idx) mc(:, idx).' * p + vc(idx, 1);
523 df_cstr = @ (p, idx, dfdp_hook) mc(:, idx).';
528 %% in a general interface, check for all(fixed) here
530 %% passed constraints
531 hook.mc = mc; % matrix of linear constraints
532 hook.vc = vc; % vector of linear constraints
533 hook.f_cstr = f_cstr; % function of all constraints
534 hook.df_cstr = df_cstr; % function of derivatives of all constraints
535 hook.n_gencstr = n_gencstr; % number of non-linear constraints
536 hook.eq_idx = false (size (vc, 1) + n_gencstr, 1);
537 hook.eq_idx(eq_idx) = true; % logical index of equality constraints in
539 hook.lbound = bounds(:, 1); % bounds, subset of linear inequality
540 % constraints in mc and vc
541 hook.ubound = bounds(:, 2);
543 %% passed values of constraints for initial parameters
544 hook.pin_cstr = pin_cstr;
546 %% passed derivative of model function
549 %% passed function for complementary pivoting
552 %% passed value of residual function for initial parameters
556 hook.max_fract_change = maxstep;
557 hook.fract_prec = pprec;
559 hook.MaxIter = niter;
561 hook.fixed = dp == 0;
563 hook.Display = 'iter';
564 __plot_cmds__ = @ __plot_cmds__; # for bug #31484 (Octave <= 3.2.4)
565 hook.plot_cmd = @ (f) __plot_cmds__ (x, y, y - f);
567 hook.Display = 'off';
570 %% only preliminary, for testing
571 hook.testing = false;
573 if (exist ('options'))
574 if (isfield (options, 'testing'))
575 hook.testing = options.testing;
577 if (isfield (options, 'new_s'))
578 hook.new_s = options.new_s;
582 [p, resid, cvg, outp] = __lm_svd__ (@ (p) y - F (x, p), pin, hook);
587 if (~cvg) disp(' CONVERGENCE NOT ACHIEVED! '); end
589 if (~(verbose || nargout > 4)) return; end
593 %% CALC VARIANCE COV MATRIX AND CORRELATION MATRIX OF PARAMETERS
594 %% re-evaluate the Jacobian at optimal values
595 jac = dFdp (p, struct ('f', f));
597 n = sum(msk); % reduce n to equal number of estimated parameters
598 jac = jac(:, msk); % use only fitted parameters
600 %% following section is Ray Muzic's estimate for covariance and correlation
601 %% assuming covariance of data is a diagonal matrix proportional to
603 %% cov matrix of data est. from Bard Eq. 7-5-13, and Row 1 Table 5.1
606 if (exist('sparse')) % save memory
607 Q = sparse (1:m, 1:m, 1 ./ tp);
608 Qinv = sparse (1:m, 1:m, tp);
610 Q = diag (ones (m, 1) ./ tp);
613 resid = resid(:); % un-weighted residuals
614 if (~isreal (resid)) error ('residuals are not real'); end
615 tp = resid.' * Qinv * resid;
616 covr = (tp / m) * Q; %covariance of residuals
618 %% Matlab compatibility and avoiding recomputation make the following
621 if (m <= n || any (eq_idx))
624 Qinv = ((m - n) / tp) * Qinv;
625 %% simplified Eq. 7-5-13, Bard; cov of parm est, inverse; outer
626 %% parantheses contain inverse of guessed covariance matrix of data
627 covpinv = jac.' * Qinv * jac;
628 if (exist ('rcond') && rcond (covpinv) <= eps)
630 elseif (rank (covpinv) < n)
631 %% above test is not equivalent to 'rcond' and may unnecessarily
632 %% reject some matrices
637 covp = inv (covpinv);
639 corp = covp ./ (d * d.');
641 covp = NA * ones (n);
646 covr=spdiags(covr,0);
648 covr=diag(covr); % convert returned values to
651 covr = reshape (covr, rows_y, cols_y);
652 stdresid = resid .* abs (wtl) / sqrt (tp / m); % equivalent to resid ./
654 stdresid = reshape (stdresid, rows_y, cols_y);
656 if (~(verbose || nargout > 8)) return; end
658 if (m > n && ~any (eq_idx))
659 Z = ((m - n) / (n * resid.' * Qinv * resid)) * covpinv;
664 %%% alt. est. of cov. mat. of parm.:(Delforge, Circulation, 82:1494-1504, 1990
665 %%disp('Alternate estimate of cov. of param. est.')
666 %%acovp=resid'*Qinv*resid/(m-n)*inv(jac'*Qinv*jac);
668 if (~(verbose || nargout > 9)) return; end
670 %%Calculate R^2, intercept form
672 tp = sumsq (yl - mean (yl));
674 r2 = 1 - sumsq (resid) / tp;
679 %% if someone has asked for it, let them have it
682 __plot_cmds__ (x, y, f);
683 disp(' Least Squares Estimates of Parameters')
685 disp(' Correlation matrix of parameters estimated')
687 disp(' Covariance matrix of Residuals' )
689 disp(' Correlation Coefficient R^2')
691 fprintf(' 95%% conf region: F(0.05)(%.0f,%.0f)>= delta_pvec.%s*Z*delta_pvec\n', n, m - n, char (39)); % works with ' and '
693 %% runs test according to Bard. p 201.
694 n1 = sum (resid > 0);
695 n2 = sum (resid < 0);
696 nrun=sum(abs(diff(resid > 0)))+1;
697 if ((n1 > 10) && (n2 > 10)) % sufficent data for test?
698 zed=(nrun-(2*n1*n2/(n1+n2)+1)+0.5)/(2*n1*n2*(2*n1*n2-n1-n2)...
699 /((n1+n2)^2*(n1+n2-1)));
701 prob = erfc(-zed/sqrt(2))/2*100;
702 disp([num2str(prob),'% chance of fewer than ',num2str(nrun),' runs.']);
704 prob = erfc(zed/sqrt(2))/2*100;
705 disp([num2str(prob),'% chance of greater than ',num2str(nrun),' runs.']);
710 function ret = tf_gencstr (p, idx, f)
712 %% necessary since user function f_gencstr might return [] or a row
718 elseif (size (ret, 2) > 1)
722 function fval = scalar_ifelse (cond, tval, fval)
724 %% needed for some anonymous functions, builtin ifelse only available
725 %% in Octave > 3.2; we need only the scalar case here
732 %! % Define functions
733 %! leasqrfunc = @(x, p) p(1) * exp (-p(2) * x);
734 %! leasqrdfdp = @(x, f, p, dp, func) [exp(-p(2)*x), -p(1)*x.*exp(-p(2)*x)];
736 %! % generate test data
739 %! data = leasqrfunc (t, p);
741 %! rnd = [0.352509; -0.040607; -1.867061; -1.561283; 1.473191; ...
742 %! 0.580767; 0.841805; 1.632203; -0.179254; 0.345208];
745 %! % wt1 = 1 /sqrt of variances of data
746 %! % 1 / wt1 = sqrt of var = standard deviation
747 %! wt1 = (1 + 0 * t) ./ sqrt (data);
748 %! data = data + 0.05 * rnd ./ wt1;
750 %! % Note by Thomas Walter <walter@pctc.chemie.uni-erlangen.de>:
752 %! % Using a step size of 1 to calculate the derivative is WRONG !!!!
753 %! % See numerical mathbooks why.
754 %! % A derivative calculated from central differences need: s
755 %! % step = 0.001...1.0e-8
756 %! % And onesided derivative needs:
757 %! % step = 1.0e-5...1.0e-8 and may be still wrong
760 %! dFdp = leasqrdfdp; % exact derivative
761 %! % dFdp = dfdp; % estimated derivative
762 %! dp = [0.001; 0.001];
764 %! stol=0.001; niter=50;
765 %! minstep = [0.01; 0.01];
766 %! maxstep = [0.8; 0.8];
767 %! options = [minstep, maxstep];
771 %! [f1, p1, kvg1, iter1, corp1, covp1, covr1, stdresid1, Z1, r21] = ...
772 %! leasqr (t, data, pin, F, stol, niter, wt1, dp, dFdp, options);
775 %! %% Example for linear inequality constraints.
776 %! %% model function:
777 %! F = @ (x, p) p(1) * exp (p(2) * x);
778 %! %% independents and dependents:
780 %! y = [1, 2, 4, 7, 14];
781 %! %% initial values:
782 %! init = [.25; .25];
783 %! %% other configuration (default values):
784 %! tolerance = .0001;
785 %! max_iterations = 20;
786 %! weights = ones (1, 5);
787 %! dp = [.001; .001]; % bidirectional numeric gradient stepsize
788 %! dFdp = 'dfdp'; % function for gradient (numerical)
790 %! %% linear constraints, A.' * parametervector + B >= 0
791 %! A = [1; -1]; B = 0; % p(1) >= p(2);
792 %! options.inequc = {A, B};
794 %! %% start leasqr, be sure that 'verbose' is not set
795 %! global verbose; verbose = false;
796 %! [f, p, cvg, iter] = ...
797 %! leasqr (x, y, init, F, tolerance, max_iterations, ...
798 %! weights, dp, dFdp, options)