Hi, I would like to make some changes (backwards-compatible, except corrections) to optim/leasqr.m (e.g. some optimization/correction in the computation of additional information like covariance matrix of parameters).
As a start, I made a diff with only style changes, to have a better base for possible discussion of later diffs. The current style is difficult to work with in Emacs' Octave-mode. I could make all changes by sending diffs to the list, but would prefer to directly change leasqr.m by svn, informing/asking the list if I think its appropriate. Regards, Olaf
Index: leasqr.m =================================================================== --- leasqr.m (revision 6632) +++ leasqr.m (working copy) @@ -1,190 +1,250 @@ -% Copyright (C) 1992-1994 Richard Shrager -% Copyright (C) 1992-1994 Arthur Jutan -% Copyright (C) 1992-1994 Ray Muzic -% -% 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 2 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/>. +## Copyright (C) 1992-1994 Richard Shrager +## Copyright (C) 1992-1994 Arthur Jutan +## Copyright (C) 1992-1994 Ray Muzic +## +## 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 2 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,kvg,iter,corp,covp,covr,stdresid,Z,r2]= ... - leasqr(x,y,pin,F,stol,niter,wt,dp,dFdp,options) -%function [f,p,kvg,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 = column vector or matrix of independent variables, 1 row per -% observation: x = [x0 x1....xm]. -% y = column vector of observed values, same number of rows as x. -% wt = column vector (dim=length(x)) of statistical weights. 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(length(y),1). -% pin = column 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; 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; 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', and 'bounds'. 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) on two -% consecutive iterations is less than their corresponding elements in -% 'options.fract_prec'. [ie. all(abs(chg*current parm est) < -% options.fract_prec) 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.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. -% -% OUTPUT VARIABLES -% f = column vector of values computed: f = F(x,p). -% p = column vector trial or final parameters. i.e, the solution. -% kvg = 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. -% -% All Zero guesses not acceptable +## -*- texinfo -*- +## +## @deftypefn {Function File} {...@var{fc}, @var{p}, @var{cvg}, @var{iter}, @var{corp}, @var{covp}, @var{covr}, @var{stdresid}, @var{Z}, @var{r2}] =} leasqr (@var{x}, @var{y}, @var{pin}, @var{F}) +## @deftypefnx {Function File} {...@var{fc}, @var{p}, @var{cvg}, @var{iter}, @var{corp}, @var{covp}, @var{covr}, @var{stdresid}, @var{Z}, @var{r2}]=} leasqr (@var{x}, @var{y}, @var{pin}, @var{F}, @var{stol}) +## @deftypefnx {Function File} {...@var{fc}, @var{p}, @var{cvg}, @var{iter}, @var{corp}, @var{covp}, @var{covr}, @var{stdresid}, @var{Z}, @var{r2}]=} leasqr (@var{x}, @var{y}, @var{pin}, @var{F}, @var{stol}, @var{niter}) +## @deftypefnx {Function File} {...@var{fc}, @var{p}, @var{cvg}, @var{iter}, @var{corp}, @var{covp}, @var{covr}, @var{stdresid}, @var{Z}, @var{r2}]=} leasqr (@var{x}, @var{y}, @var{pin}, @var{F}, @var{stol}, @var{niter}, @var{wt}) +## @deftypefnx {Function File} {...@var{fc}, @var{p}, @var{cvg}, @var{iter}, @var{corp}, @var{covp}, @var{covr}, @var{stdresid}, @var{Z}, @var{r2}]=} leasqr (@var{x}, @var{y}, @var{pin}, @var{F}, @var{stol}, @var{niter}, @var{wt}, @var{dp}) +## @deftypefnx {Function File} {...@var{fc}, @var{p}, @var{cvg}, @var{iter}, @var{corp}, @var{covp}, @var{covr}, @var{stdresid}, @var{Z}, @var{r2}]=} leasqr (@var{x}, @var{y}, @var{pin}, @var{F}, @var{stol}, @var{niter}, @var{wt}, @var{dp}, @var{dFdp}) +## @deftypefnx {Function File} {...@var{fc}, @var{p}, @var{cvg}, @var{iter}, @var{corp}, @var{covp}, @var{covr}, @var{stdresid}, @var{Z}, @var{r2}]=} leasqr (@var{x}, @var{y}, @var{pin}, @var{F}, @var{stol}, @var{niter}, @var{wt}, @var{dp}, @var{dFdp}, @var{options}) +## +## Levenberg-Marquardt nonlinear regression of @var{f}(@var{x}, @var{p}) +## to @var{y}(@var{x}). +## +## Version 3.beta +## +## @var{x} = column vector or matrix of independent variables, 1 row per +## observation: @co...@var{x} = [...@var{x0} @var{x1}....@var{xm}]}. +## +## @var{y} = column vector of observed values, same number of rows as +## @var{x}. +## +## @var{wt} = column vector (dim = length (@var{x})) of statistical +## weights. These should be set to be proportional to sqrt of @code{var +## (@var{y})^-1}; (That is, the covariance matrix of the data is assumed +## to be proportional to diagonal with diagonal equal to +## @co...@var{wt.}^-2}. The constant of proportionality will be +## estimated.); default = @code{ones (length (@var{y}), 1)}. +## +## @var{pin} = column vec of initial parameters to be adjusted by +## leasqr. +## +## @var{dp} = fractional increment of @var{p} for numerical partial +## derivatives; default = @code{.001 * ones (size (@var{pin}))}. +## @var{dp}(j) > 0 means central differences on j-th parameter +## @var{p}(j). @var{dp}(j) < 0 means one-sided differences on j-th +## parameter @var{p}(j). @var{dp}(j) = 0 holds @var{p}(j) fixed i.e. +## leasqr wont change initial guess: @var{pin}(j) +## +## @var{F} = name of function in quotes; the function shall be of the +## form +## +## @example +## +## @var{y} = f (@var{x}, @var{p}), +## +## @end example +## +## with @var{y}, @var{x}, @var{p} of the form @var{y}, @var{x}, +## @var{pin} as described above. +## +## @var{dFdp} = name of partial derivative function in quotes; default +## is "dfdp", a slow but general partial derivatives function; the +## function shall be of the form +## +## @example +## +## @var{prt} = f (@var{x}, @var{fc}, @var{p}, @var{dp}, @var{F}[, +## @var{bounds}]) . +## +## @end example +## +## For backwards compatibility, the function will only be called with an +## extra @var{bounds} argument if the @code{bounds} option is +## explicitely specified to leasqr (see dfdp.m). +## +## @var{stol} = scalar tolerance on fractional improvement in scalar sum +## of squares = @code{sum((@var{wt} .* (@var{y} - @var{f}))^2)}; default +## @var{stol} = .0001; +## +## @var{niter} = scalar maximum number of iterations; default = 20; +## +## @var{options} = structure, currently recognized fields are +## @code{fract_prec}, @code{max_fract_change}, and @code{bounds}. For +## backwards compatibility, @var{options} can also be a matrix whose +## first and second column contains the values of @code{fract_prec} and +## @code{max_fract_change}, respectively. +## +## Field @co...@var{options}.fract_prec}: column vector (same length as +## @var{pin}) of desired fractional precisions in parameter estimates. +## Iterations are terminated if change in parameter vector (@code{chg}) +## on two consecutive iterations is less than their corresponding +## elements in @var{options}.fract_prec. [ie. @code{all (abs (chg * +## current parm est) < options.fract_prec)} on two consecutive +## iterations.], default = @code{zeros ()}. +## +## Field @co...@var{options}.max_fract_change}: column vector (same +## length as @var{pin}) of maximum fractional step changes in parameter +## vector. Fractional change in elements of parameter vector is +## constrained to be at most @var{options}.max_fract_change between +## sucessive iterations. [ie. @code{abs (chg(i)) = abs (min ([chg(i) +## @var{options}.max_fract_change(i)*current param estimate])).]}, +## default = @code{Inf * ones ()}. +## +## Field @co...@var{options}.bounds}: two-column-matrix, one row for +## each parameter in @var{pin}. Each row contains a minimal and maximal +## value for each parameter. Default: @code{[-Inf, Inf]} in each row. If +## this field is used with an existing user-side function for @var{dFdp} +## (see above) the functions interface might have to be changed. +## +## OUTPUT VARIABLES +## +## @var{fc} = column vector of values computed: @co...@var{fc} = @var{F} +## (@var{x}, @var{p})}. +## +## @var{p} = column vector trial or final parameters. i.e, the solution. +## +## @var{cvg} = scalar: = 1 if convergence, = 0 otherwise. +## +## @var{iter} = scalar number of iterations used. +## +## @var{corp} = correlation matrix for parameters. +## +## @var{covp} = covariance matrix of the parameters. +## +## @var{covr} = @code{diag (covariance matrix of the residuals)}. +## +## @var{stdresid} = standardized residuals. +## +## @var{Z} = matrix that defines confidence region (see comments in the +## source). +## +## @var{r2} = coefficient of multiple determination. +## +## All Zero guesses not acceptable +## +## @end deftypefn -% 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 <[email protected]> -% Arthur Jutan <[email protected]> +function [f, p, cvg, iter, corp, covp, covr, stdresid, Z, r2] = \ + leasqr (x, y, pin, F, stol, niter,wt, dp, dFdp, options) -% 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 -% 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 bounds feature to this file and to dfdp.m, Olaf Till 4-Aug-2009 -% -% References: -% Bard, Nonlinear Parameter Estimation, Academic Press, 1974. -% Draper and Smith, Applied Regression Analysis, John Wiley and Sons, 1981. -% -%set default args + ## 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 <[email protected]> Arthur Jutan + ## <[email protected]> -% argument processing -% + ## 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 + ## 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 bounds feature to this file and to dfdp.m, Olaf Till 4-Aug-2009 + ## + ## References: + ## Bard, Nonlinear Parameter Estimation, Academic Press, 1974. + ## Draper and Smith, Applied Regression Analysis, John Wiley and Sons, 1981. + ## + ##set default args -%if (sscanf(version,'%f') >= 4), -vernum= sscanf(version,'%f'); -if vernum(1) >= 4, - global verbose - plotcmd='plot(x(:,1),y,''+'',x(:,1),f); figure(gcf)'; -else - plotcmd='plot(x(:,1),y,''+'',x(:,1),f); shg'; -end; -if (exist('OCTAVE_VERSION')) - global verbose - plotcmd='plot(x(:,1),y,"+;data;",x(:,1),f,";fit;");'; -end; + ## argument processing + ## -if(exist('verbose')~=1), %If verbose undefined, print nothing - verbose=0; %This will not tell them the results -end; + ##if (sscanf(version,'%f') >= 4), + vernum = sscanf (version, "%f"); + if (vernum(1) >= 4) + global verbose + plotcmd = 'plot (x(:,1), y, ''+'', x(:,1), f); figure (gcf)'; + else + plotcmd = 'plot (x(:,1), y, ''+'', x(:,1), f); shg'; + endif + if (exist ("OCTAVE_VERSION")) + global verbose + plotcmd= 'plot (x(:,1), y, "+;data;", x(:,1), f, ";fit;");'; + endif -if (nargin <= 8), dFdp='dfdp'; end; -if (nargin <= 7), dp=.001*(pin*0+1); end; %DT -if (nargin <= 6), wt=ones(length(y),1); end; % SMB modification -if (nargin <= 5), niter=20; end; -if (nargin == 4), stol=.0001; end; -% + if(exist ("verbose") ~= 1) # If verbose undefined, print nothing + verbose = 0; # This will not tell them the results + endif -y=y(:); wt=wt(:); pin=pin(:); dp=dp(:); %change all vectors to columns -% check data vectors- same length? -m=length(y); n=length(pin); [m1,m2]=size(x); -if m1~=m ,error('input(x)/output(y) data must have same number of rows ') ,end; + if (nargin <= 8) dFdp='dfdp'; endif + if (nargin <= 7) dp=.001*(pin*0+1); endif # DT + if (nargin <= 6) wt=ones(length(y),1); endif # SMB modification + if (nargin <= 5) niter=20; endif + if (nargin == 4) stol=.0001; endif + ## - %% processing of 'options' + y=y(:); wt=wt(:); pin=pin(:); dp=dp(:); #change all vectors to columns + ## check data vectors - same length? + m=length(y); n=length(pin); [m1,m2]=size(x); + if (m1 ~= m) + error ("input(x)/output(y) data must have same number of rows "); + endif + + ## processing of 'options' pprec = zeros (n, 1); maxstep = Inf * ones (n, 1); bounds = cat (2, -Inf * ones (n, 1), Inf * ones (n, 1)); - dfdp_cmd = "feval(dFdp,x,fbest,pprev,dp,F);"; % will possibly be redefined + dfdp_cmd = "feval(dFdp,x,fbest,pprev,dp,F);"; # will possibly be redefined if (nargin > 9) - if (ismatrix (options)) % backwards compatibility + if (ismatrix (options)) # backwards compatibility tp = options; options = struct ("fract_prec", tp(:, 1)); if (columns (tp) > 1) @@ -232,195 +292,201 @@ error ("no free parameters"); endif -% set up for iterations -% -p = pin; -f=feval(F,x,p); fbest=f; pbest=p; -r=wt.*(y-f); -ss=r'*r; -sbest=ss; -nrm=zeros(n,1); -chgprev=Inf*ones(n,1); -kvg=0; -epsLlast=1; -epstab=[.1, 1, 1e2, 1e4, 1e6]; + ## set up for iterations + ## + p = pin; + f = feval (F, x, p); fbest = f; pbest = p; + r = wt .* (y - f); + ss = r' * r; + sbest = ss; + nrm = zeros (n, 1); + chgprev = Inf * ones (n, 1); + cvg = 0; + epsLlast = 1; + epstab = [.1, 1, 1e2, 1e4, 1e6]; -% do iterations -% -for iter=1:niter, - pprev=pbest; - prt = eval (dfdp_cmd); - r=wt.*(y-fbest); - sprev=sbest; - sgoal=(1-stol)*sprev; - for j=1:n, - if dp(j)==0, - nrm(j)=0; + ## do iterations + ## + for iter = 1:niter + pprev = pbest; + prt = eval (dfdp_cmd); + r = wt .* (y - fbest); + sprev = sbest; + sgoal = (1 - stol) * sprev; + for j=1:n + if (dp (j) == 0) + nrm(j) = 0; + else + prt (:, j) = wt .* prt (:, j); + nrm (j) = prt(:, j)' * prt(:,j); + if (nrm(j) > 0) + nrm(j) = 1 / sqrt (nrm(j)); + endif + endif + prt(:,j) = nrm(j) * prt(:,j); + endfor + ## above loop could ? be replaced by: prt=prt.*wt(:,ones(1,n)); + ## nrm=dp./sqrt(diag(prt'*prt)); prt=prt.*nrm(:,ones(1,m))'; + [prt, s, v] = svd (prt,0); + s = diag (s); + g = prt' * r; + for jjj = 1:length (epstab) + epsL = max (epsLlast * epstab(jjj), 1e-7); + se = sqrt ((s .* s) + epsL); + gse = g ./ se; + chg = ((v * gse) .* nrm); + ## check the change constraints and apply as necessary + ochg = chg; + idx = ~isinf (maxstep); + limit = abs (maxstep(idx) .* pprev(idx)); + chg(idx) = min(max(chg(idx), -limit), limit); + if (verbose && any(ochg ~= chg)) + disp (["Change in parameter(s): ", \ + sprintf("%d ", find(ochg ~= chg)), "were constrained"]); + endif + aprec = abs (pprec .* pbest); # --- + ## ss=scalar sum of squares=sum((wt.*(y-f))^2). + if (any (abs (chg) > 0.1 * aprec)) # --- # only worth evaluating + # function if there is some + # non-miniscule change + p = chg + pprev; + ## apply bounds; preserving the direction of the attempted step + ## might lead to fixing _all_ parameters at their current value, + ## so decided not to do that and to simply reset parameters to + ## bounds + lidx = p < bounds(:, 1); + p(lidx) = bounds(lidx, 1); + uidx = p > bounds(:, 2); + p(uidx) = bounds(uidx, 2); + if (verbose && any (idx = lidx | uidx)) + printf ("bounds apply for parameters number %s %i\n", \ + sprintf ("%i, ", find (idx)(1:end - 1)), \ + find (idx)(end)); + endif + ## + f = feval ( F, x, p); + r = wt .* (y - f); + ss = r' * r; + if (ss < sbest) + pbest = p; + fbest = f; + sbest = ss; + endif + if (ss <= sgoal) + break; + endif + endif #--- + endfor + epsLlast = epsL; + if (verbose) + eval (plotcmd); + endif + if (ss < eps) + break; + endif + aprec = abs (pprec .* pbest); + ## [aprec, chg, chgprev] + if (all (abs (chg) < aprec) && all (abs (chgprev) < aprec)) + cvg = 1; + if (verbose) + fprintf ("Parameter changes converged to specified precision\n"); + endif + break; else - prt(:,j)=wt.*prt(:,j); - nrm(j)=prt(:,j)'*prt(:,j); - if nrm(j)>0, - nrm(j)=1/sqrt(nrm(j)); - end; - end - prt(:,j)=nrm(j)*prt(:,j); - end; -% above loop could ? be replaced by: -% prt=prt.*wt(:,ones(1,n)); -% nrm=dp./sqrt(diag(prt'*prt)); -% prt=prt.*nrm(:,ones(1,m))'; - [prt,s,v]=svd(prt,0); - s=diag(s); - g=prt'*r; - for jjj=1:length(epstab), - epsL = max(epsLlast*epstab(jjj),1e-7); - se=sqrt((s.*s)+epsL); - gse=g./se; - chg=((v*gse).*nrm); -% check the change constraints and apply as necessary - ochg=chg; - idx = ~isinf(maxstep); - limit = abs(maxstep(idx).*pprev(idx)); - chg(idx) = min(max(chg(idx),-limit),limit); - if (verbose & any(ochg ~= chg)), - disp(['Change in parameter(s): ', ... - sprintf('%d ',find(ochg ~= chg)), 'were constrained']); - end; - aprec=abs(pprec.*pbest); %--- -% ss=scalar sum of squares=sum((wt.*(y-f))^2). - if (any(abs(chg) > 0.1*aprec)),%--- % only worth evaluating function if - p=chg+pprev; % there is some non-miniscule change - %% apply bounds; preserving the direction of the attempted step - %% might lead to fixing _all_ parameters at their current value, - %% so decided not to do that and to simply reset parameters to - %% bounds - lidx = p < bounds(:, 1); - p(lidx) = bounds(lidx, 1); - uidx = p > bounds(:, 2); - p(uidx) = bounds(uidx, 2); - if (verbose && any (idx = lidx | uidx)) - printf ("bounds apply for parameters number %s %i\n", \ - sprintf ("%i, ", find (idx)(1:end - 1)), \ - find (idx)(end)); - endif - %% - f=feval(F,x,p); - r=wt.*(y-f); - ss=r'*r; - if ss<sbest, - pbest=p; - fbest=f; - sbest=ss; - end; - if ss<=sgoal, - break; - end; - end; %--- - end; - epsLlast = epsL; - if (verbose), - eval(plotcmd); - end; - if ss<eps, - break; - end - aprec=abs(pprec.*pbest); -% [aprec, chg, chgprev] - if (all(abs(chg) < aprec) & all(abs(chgprev) < aprec)), - kvg=1; - if (verbose), - fprintf('Parameter changes converged to specified precision\n'); - end; - break; - else - chgprev=chg; - end; - if ss>sgoal, - break; - end; -end; + chgprev = chg; + endif + if (ss > sgoal) + break; + endif + endfor -% set return values -% -p=pbest; -f=fbest; -ss=sbest; -kvg=((sbest>sgoal)|(sbest<=eps)|kvg); -if kvg ~= 1 , disp(' CONVERGENCE NOT ACHIEVED! '), end; + ## set return values + ## + p = pbest; + f = fbest; + ss = sbest; + cvg = ((sbest > sgoal) || (sbest <= eps) || cvg); + if (cvg ~= 1) disp (" CONVERGENCE NOT ACHIEVED! "); endif -% CALC VARIANCE COV MATRIX AND CORRELATION MATRIX OF PARAMETERS -% re-evaluate the Jacobian at optimal values -jac = eval (dfdp_cmd); -msk = dp ~= 0; -n = sum(msk); % reduce n to equal number of estimated parameters -jac = jac(:, msk); % use only fitted parameters + ## CALC VARIANCE COV MATRIX AND CORRELATION MATRIX OF PARAMETERS + ## re-evaluate the Jacobian at optimal values + jac = eval (dfdp_cmd); + msk = dp ~= 0; + 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 + ## 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 -if exist('sparse') % save memory - Q=sparse(1:m,1:m,1./wt.^2); - Qinv=sparse(1:m,1:m,wt.^2); -else - Q=diag((0*wt+1)./(wt.^2)); - Qinv=diag(wt.*wt); -end -resid=y-f; %un-weighted residuals -covr=resid'*Qinv*resid*Q/(m-n); %covariance of residuals -Vy=1/(1-n/m)*covr; % Eq. 7-13-22, Bard %covariance of the data + if (exist ("sparse")) # save memory + Q = sparse (1:m, 1:m, 1 ./ wt.^2); + Qinv = sparse (1:m, 1:m, wt.^2); + else + Q = diag ((0 * wt + 1) ./ (wt.^2)); + Qinv = diag (wt .* wt); + endif + resid = y - f; # un-weighted residuals + covr = resid' * Qinv * resid * Q / (m - n); # covariance of residuals + Vy = 1 / (1 - n / m) * covr; # Eq. 7-13-22, Bard # covariance of the data -jtgjinv=inv(jac'*Qinv*jac); %argument of inv may be singular -covp=jtgjinv*jac'*Qinv*Vy*Qinv*jac*jtgjinv; % Eq. 7-5-13, Bard %cov of parm est -d=sqrt(abs(diag(covp))); -corp=covp./(d*d'); + jtgjinv = inv (jac' * Qinv * jac); # argument of inv may be singular + covp = jtgjinv * jac' * Qinv * Vy * Qinv * jac * jtgjinv; #Eq.7-5-13, + #Bard #cov of parm est + d = sqrt (abs (diag (covp))); + corp = covp ./ (d * d'); -if exist('sparse') - covr=spdiags(covr,0); - stdresid=resid./sqrt(spdiags(Vy,0)); -else - covr=diag(covr); % convert returned values to compact storage - stdresid=resid./sqrt(diag(Vy)); % compute then convert for compact storage -end -Z=((m-n)*jac'*Qinv*jac)/(n*resid'*Qinv*resid); + if (exist("sparse")) + covr=spdiags(covr,0); + stdresid=resid./sqrt(spdiags(Vy,0)); + else + covr=diag(covr); # convert returned values to compact storage + stdresid=resid./sqrt(diag(Vy)); # compute then convert for compact storage + endif + Z = ((m - n) * jac' * Qinv * jac) / (n * resid' * Qinv * resid); -%%% 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)*jtgjinv +### 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)*jtgjinv -%Calculate R^2 (Ref Draper & Smith p.46) -% -r=corrcoef([y(:),f(:)]); -r2=r(1,2).^2; + ##Calculate R^2 (Ref Draper & Smith p.46) + ## + r = corrcoef ([y(:), f(:)]); + r2 = r(1, 2).^2; -% if someone has asked for it, let them have it -% -if (verbose), - eval(plotcmd); - 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) - sprintf(' 95%% conf region: F(0.05)(%.0f,%.0f)>= delta_pvec''*Z*delta_pvec',n,m-n) - Z -% runs test according to Bard. p 201. - n1 = sum((f-y) < 0); - n2 = sum((f-y) > 0); - nrun=sum(abs(diff((f-y)<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; + ## if someone has asked for it, let them have it + ## + if (verbose), + eval (plotcmd); + 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) + sprintf (" 95%% conf region: F(0.05)(%.0f,%.0f) >= delta_pvec' * Z * delta_pvec", \ + n, m - n) + Z + ## runs test according to Bard. p 201. + n1 = sum ((f - y) < 0); + n2 = sum ((f - y) > 0); + nrun = sum (abs (diff ((f - y) < 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."]); + endif + endif + endif
------------------------------------------------------------------------------ Return on Information: Google Enterprise Search pays you back Get the facts. http://p.sf.net/sfu/google-dev2dev
_______________________________________________ Octave-dev mailing list [email protected] https://lists.sourceforge.net/lists/listinfo/octave-dev
