On Tue, Aug 04, 2009 at 09:57:28PM +0200, Søren Hauberg wrote:
> tir, 04 08 2009 kl. 13:02 +0200, skrev Olaf Till:
> > A lot of
> > cosmetic changes have been made to make leasqr.m better editable under
> > emacs, so almost each line is changed and its not much use making
> > diffs :-(.
>
> Then could you please start out by doing a diff where you do not change
> coding style. Once we have the content in place, you can always change
> the coding style.
Ok, attached is an SVN diff inside main/optim/inst/.
Olaf
Index: leasqr.m
===================================================================
--- leasqr.m (revision 6076)
+++ leasqr.m (working copy)
@@ -46,17 +46,28 @@
% 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 = matrix of n rows (same number of rows as pin) containing
-% column 1: desired fractional precision in parameter estimates.
-% Iterations are terminated if change in parameter vector (chg) on two
-% consecutive iterations is less than their corresponding elements
-% in options(:,1). [ie. all(abs(chg*current parm est) < options(:,1))
-% on two consecutive iterations.], default = zeros().
-% column 2: maximum fractional step change in parameter vector.
-% Fractional change in elements of parameter vector is constrained to be
-% at most options(:,2) between sucessive iterations.
-% [ie. abs(chg(i))=abs(min([chg(i) options(i,2)*current param estimate])).],
-% default = Inf*ones().
+% 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.
%
% OUTPUT VARIABLES
% f = column vector of values computed: f = F(x,p).
@@ -122,6 +133,7 @@
% 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.
@@ -158,27 +170,64 @@
y=y(:); wt=wt(:); pin=pin(:); dp=dp(:); %change all vectors to columns
% check data vectors- same length?
-m=length(y); n=length(pin); p=pin;[m1,m2]=size(x);
+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 <= 9),
- options=[zeros(n,1), Inf*ones(n,1)];
- nor = n; noc = 2;
-else
- [nor, noc]=size(options);
- if (nor ~= n),
- error('options and parameter matrices must have same number of rows'),
- end;
- if (noc ~= 2),
- options=[options(:,1), Inf*ones(nor,1)];
- end;
-end;
-pprec=options(:,1);
-maxstep=options(:,2);
-%
+ %% processing of 'options'
+ pprec = zeros (n, 1);
+ maxstep = Inf * ones (n, 1);
+ bounds = cat (2, -Inf * ones (n, 1), Inf * ones (n, 1));
+ 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);
+ endif
+ endif
+ if (isfield (options, "fract_prec"))
+ pprec = options.fract_prec;
+ if (rows (pprec) != n || columns (pprec) != 1)
+ error ("fractional precisions: wrong dimensions");
+ endif
+ endif
+ if (isfield (options, "max_fract_change"))
+ maxstep = options.max_fract_change;
+ if (rows (maxstep) != n || columns (maxstep) != 1)
+ error ("maximum fractional step changes: wrong dimensions");
+ endif
+ endif
+ if (isfield (options, "bounds"))
+ bounds = options.bounds;
+ if (rows (bounds) != n || columns (bounds) != 2)
+ error ("bounds: wrong dimensions");
+ endif
+ idx = bounds(:, 1) > bounds(:, 2);
+ tp = bounds(idx, 2);
+ bounds(idx, 2) = bounds(idx, 1);
+ bounds(idx, 1) = tp;
+ if (any (idx = bounds(:, 1) == bounds(:, 2)))
+ warning ("lower and upper bounds identical for some parameters, setting the respective elements of 'dp' to zero");
+ dp(idx) = 0;
+ endif
+ if (any (idx = pin < bounds(:, 1)))
+ warning ("some initial parameters set to lower bound");
+ pin(idx) = bounds(idx, 1);
+ endif
+ if (any (idx = pin > bounds(:, 2)))
+ warning ("some initial parameters set to upper bound");
+ pin(idx) = bounds(idx, 2);
+ endif
+ endif
+ endif
+ if (all (dp == 0))
+ 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;
@@ -193,7 +242,7 @@
%
for iter=1:niter,
pprev=pbest;
- prt=feval(dFdp,x,fbest,pprev,dp,F);
+ prt=feval(dFdp,x,fbest,pprev,dp,F,bounds);
r=wt.*(y-fbest);
sprev=sbest;
sgoal=(1-stol)*sprev;
@@ -234,6 +283,20 @@
% 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;
@@ -280,7 +343,7 @@
% CALC VARIANCE COV MATRIX AND CORRELATION MATRIX OF PARAMETERS
% re-evaluate the Jacobian at optimal values
-jac=feval(dFdp,x,f,p,dp,F);
+jac=feval(dFdp,x,f,p,dp,F,bounds);
msk = dp ~= 0;
n = sum(msk); % reduce n to equal number of estimated parameters
jac = jac(:, msk); % use only fitted parameters
Index: dfdp.m
===================================================================
--- dfdp.m (revision 6076)
+++ dfdp.m (working copy)
@@ -15,7 +15,7 @@
% 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 prt=dfdp(x,f,p,dp,func)
+function prt=dfdp(x,f,p,dp,func,bounds)
% numerical partial derivatives (Jacobian) df/dp for use with leasqr
% --------INPUT VARIABLES---------
% x=vec or matrix of indep var(used as arg to func) x=[x0 x1 ....]
@@ -27,26 +27,55 @@
% dp(j)=0 sets corresponding partials to zero; i.e. holds p(j) fixed
% func=string naming the function (.m) file
% e.g. to calc Jacobian for function expsum prt=dfdp(x,f,p,dp,'expsum')
+% bounds=two-column-matrix of lower and upper bounds for parameters
%----------OUTPUT VARIABLES-------
% prt= Jacobian Matrix prt(i,j)=df(i)/dp(j)
%================================
m=size(x,1); if (m==1), m=size(x,2); end %# PAK: in case #cols > #rows
n=length(p); %dimensions
-ps=p; prt=zeros(m,n);del=zeros(n,1); % initialise Jacobian to Zero
-for j=1:n
- del(j)=dp(j) .*p(j); %cal delx=fract(dp)*param value(p)
- if p(j)==0
- del(j)=dp(j); %if param=0 delx=fraction
- end
- p(j)=ps(j) + del(j);
- if del(j)~=0, f1=feval(func,x,p);
- if dp(j) < 0, prt(:,j)=(f1-f)./del(j);
- else
- p(j)=ps(j)- del(j);
- prt(:,j)=(f1-feval(func,x,p))./(2 .*del(j));
- end
- end
- p(j)=ps(j); %restore p(j)
-end
-return
+prt=zeros(m,n);del=zeros(n,1); % initialise Jacobian to Zero
+del = dp .* p; %cal delx=fract(dp)*param value(p)
+idx = p == 0;
+del(idx) = dp(idx); %if param=0 delx=fraction
+idx = dp > 0;
+del(idx) = abs (del(idx)); % not for one-sided intervals, changed
+ % direction of intervals could change
+ % behavior of optimization without bounds
+min_del = min (abs (del), bounds(:, 2) - bounds(:, 1));
+ for j=1:n
+ ps = p;
+ if (dp(j)~=0)
+ if (dp(j) < 0)
+ ps(j) = p(j) + del(j);
+ if (ps(j) < bounds(j, 1) || ps(j) > bounds(j, 2))
+ t_del1 = max (bounds(j, 1) - p(j), - abs (del(j))); %
+ %non-positive
+ t_del2 = min (bounds(j, 2) - p(j), abs (del(j))); %
+ %non-negative
+ if (- t_del1 > t_del2)
+ del(j) = t_del1;
+ else
+ del(j) = t_del2;
+ endif
+ ps(j) = p(j) + del(j);
+ endif
+ prt(:, j) = (feval (func, x, ps) - f) / del(j);
+ else
+ if (p(j) - del(j) < bounds(j, 1))
+ tp = bounds(j, 1);
+ ps(j) = tp + min_del(j);
+ elseif (p(j) + del(j) > bounds(j, 2))
+ ps(j) = bounds(j, 2);
+ tp = ps(j) - min_del(j);
+ else
+ ps(j) = p(j) + del(j);
+ tp = p(j) - del(j);
+ min_del(j) = 2 * del(j);
+ endif
+ f1 = feval (func, x, ps);
+ ps(j) = tp;
+ prt(:, j) = (f1 - feval (func, x, ps)) / (min_del(j));
+ endif
+ endif
+ endfor
------------------------------------------------------------------------------
Let Crystal Reports handle the reporting - Free Crystal Reports 2008 30-Day
trial. Simplify your report design, integration and deployment - and focus on
what you do best, core application coding. Discover what's new with
Crystal Reports now. http://p.sf.net/sfu/bobj-july
_______________________________________________
Octave-dev mailing list
Octave-dev@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/octave-dev
