Mercurial > forge

comparison main/optim/inst/leasqr.m @ 9930:d30cfca46e8a octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
optim: upgrade license to GPLv3+ and mention on DESCRIPTION the other package licenses
author carandraug
date Fri, 30 Mar 2012 15:14:48 +0000
parents 88dc7961d556
children fba8cdd5f9ad
comparison
equal deleted inserted replaced
9929:df50d0ae107f 9930:d30cfca46e8a
1 %% Copyright (C) 1992-1994 Richard Shrager 1 %% Copyright (C) 1992-1994 Richard Shrager
2 %% Copyright (C) 1992-1994 Arthur Jutan 2 %% Copyright (C) 1992-1994 Arthur Jutan <jutan@charon.engga.uwo.ca>
3 %% Copyright (C) 1992-1994 Ray Muzic 3 %% Copyright (C) 1992-1994 Ray Muzic <rfm2@ds2.uh.cwru.edu>
4 %% Copyright (C) 2010, 2011 Olaf Till <olaf.till@uni-jena.de> 4 %% Copyright (C) 2010, 2011 Olaf Till <olaf.till@uni-jena.de>
5 %% 5 %%
6 %% This program is free software; you can redistribute it and/or modify 6 %% This program is free software; you can redistribute it and/or modify it under
7 %% it under the terms of the GNU General Public License as published by 7 %% the terms of the GNU General Public License as published by the Free Software
8 %% the Free Software Foundation; either version 2 of the License, or 8 %% Foundation; either version 3 of the License, or (at your option) any later
9 %% (at your option) any later version. 9 %% version.
10 %% 10 %%
11 %% This program is distributed in the hope that it will be useful, 11 %% This program is distributed in the hope that it will be useful, but WITHOUT
12 %% but WITHOUT ANY WARRANTY; without even the implied warranty of 12 %% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 %% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 %% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
14 %% GNU General Public License for more details. 14 %% details.
15 %% 15 %%
16 %% You should have received a copy of the GNU General Public License 16 %% You should have received a copy of the GNU General Public License along with
17 %% along with this program; If not, see <http://www.gnu.org/licenses/>. 17 %% this program; if not, see <http://www.gnu.org/licenses/>.
18
19 %%function [f,p,cvg,iter,corp,covp,covr,stdresid,Z,r2]=
20 %% leasqr(x,y,pin,F,{stol,niter,wt,dp,dFdp,options})
21 %%
22 %% Levenberg-Marquardt nonlinear regression of f(x,p) to y(x).
23 %%
24 %% Version 3.beta
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 =
69 %% Inf*ones().
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.
132 %%
133 %% OUTPUT VARIABLES
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.
144 %%
145 %% Not suitable for non-real residuals.
146 %%
147 %% References:
148 %% Bard, Nonlinear Parameter Estimation, Academic Press, 1974.
149 %% Draper and Smith, Applied Regression Analysis, John Wiley and Sons, 1981.
18 150
19 function [f,p,cvg,iter,corp,covp,covr,stdresid,Z,r2]= ... 151 function [f,p,cvg,iter,corp,covp,covr,stdresid,Z,r2]= ...
20 leasqr(x,y,pin,F,stol,niter,wt,dp,dFdp,options) 152 leasqr(x,y,pin,F,stol,niter,wt,dp,dFdp,options)
21
22 %%function [f,p,cvg,iter,corp,covp,covr,stdresid,Z,r2]=
23 %% leasqr(x,y,pin,F,{stol,niter,wt,dp,dFdp,options})
24 %%
25 %% Levenberg-Marquardt nonlinear regression of f(x,p) to y(x).
26 %%
27 %% Version 3.beta
28 %% Optional parameters are in braces {}.
29 %% x = vector or matrix of independent variables.
30 %% y = vector or matrix of observed values.
31 %% wt = statistical weights (same dimensions as y). These should be
32 %% set to be proportional to (sqrt of var(y))^-1; (That is, the
33 %% covariance matrix of the data is assumed to be proportional to
34 %% diagonal with diagonal equal to (wt.^2)^-1. The constant of
35 %% proportionality will be estimated.); default = ones( size (y)).
36 %% pin = vec of initial parameters to be adjusted by leasqr.
37 %% dp = fractional increment of p for numerical partial derivatives;
38 %% default = .001*ones(size(pin))
39 %% dp(j) > 0 means central differences on j-th parameter p(j).
40 %% dp(j) < 0 means one-sided differences on j-th parameter p(j).
41 %% dp(j) = 0 holds p(j) fixed i.e. leasqr wont change initial guess: pin(j)
42 %% F = name of function in quotes or function handle; the function
43 %% shall be of the form y=f(x,p), with y, x, p of the form y, x, pin
44 %% as described above.
45 %% dFdp = name of partial derivative function in quotes or function
46 %% handle; default is 'dfdp', a slow but general partial derivatives
47 %% function; the function shall be of the form
48 %% prt=dfdp(x,f,p,dp,F[,bounds]). For backwards compatibility, the
49 %% function will only be called with an extra 'bounds' argument if the
50 %% 'bounds' option is explicitely specified to leasqr (see dfdp.m).
51 %% stol = scalar tolerance on fractional improvement in scalar sum of
52 %% squares = sum((wt.*(y-f))^2); default stol = .0001;
53 %% niter = scalar maximum number of iterations; default = 20;
54 %% options = structure, currently recognized fields are 'fract_prec',
55 %% 'max_fract_change', 'inequc', 'bounds', and 'equc'. For backwards
56 %% compatibility, 'options' can also be a matrix whose first and
57 %% second column contains the values of 'fract_prec' and
58 %% 'max_fract_change', respectively.
59 %% Field 'options.fract_prec': column vector (same length as 'pin')
60 %% of desired fractional precisions in parameter estimates.
61 %% Iterations are terminated if change in parameter vector (chg)
62 %% relative to current parameter estimate is less than their
63 %% corresponding elements in 'options.fract_prec' [ie. all (abs
64 %% (chg) < abs (options.fract_prec .* current_parm_est))] on two
65 %% consecutive iterations, default = zeros().
66 %% Field 'options.max_fract_change': column vector (same length as
67 %% 'pin) of maximum fractional step changes in parameter vector.
68 %% Fractional change in elements of parameter vector is constrained to
69 %% be at most 'options.max_fract_change' between sucessive iterations.
70 %% [ie. abs(chg(i))=abs(min([chg(i)
71 %% options.max_fract_change(i)*current param estimate])).], default =
72 %% Inf*ones().
73 %% Field 'options.inequc': cell-array containing up to four entries,
74 %% two entries for linear inequality constraints and/or one or two
75 %% entries for general inequality constraints. Initial parameters
76 %% must satisfy these constraints. Either linear or general
77 %% constraints may be the first entries, but the two entries for
78 %% linear constraints must be adjacent and, if two entries are given
79 %% for general constraints, they also must be adjacent. The two
80 %% entries for linear constraints are a matrix (say m) and a vector
81 %% (say v), specifying linear inequality constraints of the form
82 %% `m.' * parameters + v >= 0'. If the constraints are just bounds,
83 %% it is suggested to specify them in 'options.bounds' instead,
84 %% since then some sanity tests are performed, and since the
85 %% function 'dfdp.m' is guarantied not to violate constraints during
86 %% determination of the numeric gradient only for those constraints
87 %% specified as 'bounds' (possibly with violations due to a certain
88 %% inaccuracy, however, except if no constraints except bounds are
89 %% specified). The first entry for general constraints must be a
90 %% differentiable vector valued function (say h), specifying general
91 %% inequality constraints of the form `h (p[, idx]) >= 0'; p is the
92 %% column vector of optimized paraters and the optional argument idx
93 %% is a logical index. h has to return the values of all constraints
94 %% if idx is not given, and has to return only the indexed
95 %% constraints if idx is given (so computation of the other
96 %% constraints can be spared). If a second entry for general
97 %% constraints is given, it must be a function (say dh) which
98 %% returnes a matrix whos rows contain the gradients of the
99 %% constraint function h with respect to the optimized parameters.
100 %% It has the form jac_h = dh (vh, p, dp, h, idx[, bounds]); p is
101 %% the column vector of optimized parameters, and idx is a logical
102 %% index --- only the rows indexed by idx must be returned (so
103 %% computation of the others can be spared). The other arguments of
104 %% dh are for the case that dh computes numerical gradients: vh is
105 %% the column vector of the current values of the constraint
106 %% function h, with idx already applied. h is a function h (p) to
107 %% compute the values of the constraints for parameters p, it will
108 %% return only the values indexed by idx. dp is a suggestion for
109 %% relative step width, having the same value as the argument 'dp'
110 %% of leasqr above. If bounds were specified to leasqr, they are
111 %% provided in the argument bounds of dh, to enable their
112 %% consideration in determination of numerical gradients. If dh is
113 %% not specified to leasqr, numerical gradients are computed in the
114 %% same way as with 'dfdp.m' (see above). If some constraints are
115 %% linear, they should be specified as linear constraints (or
116 %% bounds, if applicable) for reasons of performance, even if
117 %% general constraints are also specified.
118 %% Field 'options.bounds': two-column-matrix, one row for each
119 %% parameter in 'pin'. Each row contains a minimal and maximal value
120 %% for each parameter. Default: [-Inf, Inf] in each row. If this
121 %% field is used with an existing user-side function for 'dFdp'
122 %% (see above) the functions interface might have to be changed.
123 %% Field 'options.equc': equality constraints, specified the same
124 %% way as inequality constraints (see field 'options.inequc').
125 %% Initial parameters must satisfy these constraints.
126 %% Note that there is possibly a certain inaccuracy in honoring
127 %% constraints, except if only bounds are specified.
128 %% _Warning_: If constraints (or bounds) are set, returned guesses
129 %% of corp, covp, and Z are generally invalid, even if no constraints
130 %% are active for the final parameters. If equality constraints are
131 %% specified, corp, covp, and Z are not guessed at all.
132 %% Field 'options.cpiv': Function for complementary pivot algorithm
133 %% for inequality constraints, default: cpiv_bard. No different
134 %% function is supplied.
135 %%
136 %% OUTPUT VARIABLES
137 %% f = column vector of values computed: f = F(x,p).
138 %% p = column vector trial or final parameters. i.e, the solution.
139 %% cvg = scalar: = 1 if convergence, = 0 otherwise.
140 %% iter = scalar number of iterations used.
141 %% corp = correlation matrix for parameters.
142 %% covp = covariance matrix of the parameters.
143 %% covr = diag(covariance matrix of the residuals).
144 %% stdresid = standardized residuals.
145 %% Z = matrix that defines confidence region (see comments in the source).
146 %% r2 = coefficient of multiple determination, intercept form.
147 %%
148 %% Not suitable for non-real residuals.
149 153
150 %% The following two blocks of comments are chiefly from the original 154 %% The following two blocks of comments are chiefly from the original
151 %% version for Matlab. For later changes the logs of the Octave Forge 155 %% version for Matlab. For later changes the logs of the Octave Forge
152 %% svn repository should also be consulted. 156 %% svn repository should also be consulted.
153 157
215 %% overall cycle with a new epsL, some performance gains from linear 219 %% overall cycle with a new epsL, some performance gains from linear
216 %% constraints even if general constraints are specified. Equality 220 %% constraints even if general constraints are specified. Equality
217 %% constraints also implemented. Olaf Till 221 %% constraints also implemented. Olaf Till
218 %% Now split into files leasqr.m and __lm_svd__.m. 222 %% Now split into files leasqr.m and __lm_svd__.m.
219 223
220 %%
221 %% References:
222 %% Bard, Nonlinear Parameter Estimation, Academic Press, 1974.
223 %% Draper and Smith, Applied Regression Analysis, John Wiley and Sons, 1981.
224
225 %% needed for some anonymous functions 224 %% needed for some anonymous functions
226 if (exist ('ifelse') ~= 5) 225 if (exist ('ifelse') ~= 5)
227 ifelse = @ scalar_ifelse; 226 ifelse = @ scalar_ifelse;
228 end 227 end
229 228