Hi All,
I have created and modified two algorithms.
One of them is the nonlinear conjugate gradient explained on the
wiki page. The code is based on Etienne Grossmann code however
I had to modify it seriously.
Now it is in sync with the wikipedia page description.
There are two comments I have to make.
Etienne said that BGFS method works better in many cases
which for this reason is already called in
fmins and fsearchmin functions.
These algorithms do not rely on analytical derivatives
which may or many not be an advantage.
The cg_min function I provided requires analytical derivatives and the main
purpose to provide it is to have it since many articles do performance studies
with
nonlin conjgrad.
The other algo of mine is written from scratch.
It gives the confidence level of being a multinomial distribution well
represented by
the sampling portions (ie. within predefined confidence intervals ).
It considers finite population and infinite populations as well.
Please give comments on my functions and let me put in the repository.
Levente
optim package:
cg_min.m
bracket_min.m
brent_line_min.m
semi_bracket.m
statistics package:
cl_multinom.m
Thanks:
Levente
--
Blogger of http://fapuma.blogspot.com
Pictures of my life at http://picasaweb.google.com/TorokLev
In addition: http://www.youtube.com/watch?v=0CUDg3NPEXY&fmt=16
-------------------------------------------------------
--
Blogger of http://fapuma.blogspot.com
## -*- texinfo -*-
## @deftypefn {Function File} {...@var{s},@var{v},@var{n}]} brent_line_min ( @var{f},@var{df},@var{args},@var{ctl} )
## Line minimization of f along df
##
## Finds minimum of f on line @math{ x0 + dx*w | a < w < b } by
## bracketing. a and b are passed through argument ctl.
##
## @subheading Arguments
## @itemize @bullet
## @item @var{f} : string : Name of function. Must return a real value
## @item @var{args} : cell : Arguments passed to f or RxC : f's only argument. x0 must be at @var{ar...@{ @var{ctl}(2) @}
## @item @var{ctl} : 5 : (optional) Control variables, described below.
## @end itemize
##
## @subheading Returned values
## @itemize @bullet
## @item @var{s} : 1 : Minimum is at x0 + s*dx
## @item @var{v} : 1 : Value of f at x0 + s*dx
## @item @var{nev} : 1 : Number of function evaluations
## @end itemize
##
## @subheading Control Variables
## @itemize @bullet
## @item @var{ctl}(1) : Upper bound for error on s Default=sqrt(eps)
## @item @var{ctl}(2) : Position of minimized argument in args Default= 1
## @item @var{ctl}(3) : Maximum number of function evaluations Default= inf
## @item @var{ctl}(4) : a Default=-inf
## @item @var{ctl}(5) : b Default= inf
## @end itemize
##
## Default values will be used if ctl is not passed or if nan values are
## given.
## @end deftypefn
function [s,gs,nev] = brent_line_min( f,dx,args,ctl )
verbose = 0;
seps = sqrt (eps);
# Default control variables
tol = 10*eps; # sqrt (eps);
narg = 1;
maxev = inf;
a = -inf;
b = inf;
if nargin >= 4, # Read arguments
if !isnan (ctl (1)), tol = ctl(1); end
if length (ctl)>=2 && !isnan (ctl(2)), narg = ctl(2); end
if length (ctl)>=3 && !isnan (ctl(3)), maxev = ctl(3); end
if length (ctl)>=4 && !isnan (ctl(4)), a = ctl(4); end
if length (ctl)>=5 && !isnan (ctl(5)), b = ctl(5); end
end # Otherwise, use defaults, def'd above
if a>b, tmp=a; a=b; b=tmp; end
if narg > length (args),
printf ("brent_line_min : narg==%i > length (args)==%i",\
narg, length (args));
keyboard
end
if ! iscell (args),
args = {args};
endif
x = args{ narg };
[R,C] = size (x);
N = R*C; # Size of argument
gs0 = gs = feval (f, args);
nev = 1;
# Initial value
s = 0;
if all (dx==0), return; end # Trivial case
# If any of the bounds is infinite, find
# finite bounds that bracket minimum
if !isfinite (a) || !isfinite (b),
if !isfinite (a) && !isfinite (b),
[a,b, ga,gb, n] = bracket_min (f, dx, narg, args);
elseif !isfinite (a),
[a,b, ga,gb, n] = semi_bracket (f, -dx, -b, narg, args);
tmp = a; a = -b; b = -tmp;
tmp = ga; ga = gb; gb = tmp;
else
[a,b, ga,gb, n] = semi_bracket (f, dx, a, narg, args);
end
nev += n;
else
args{narg} = x+a*dx; ga = feval( f, args );
args{narg} = x+b*dx; gb = feval( f, args );
nev += 2;
end # End of finding bracket for minimum
if a > b, # Check assumptions
printf ("brent_line_min : a > b\n");
keyboard
end
s = 0.5*(a+b);
args{narg} = x+ s*dx; gs = feval( f, args );
end
nev++;
if verbose,
printf ("[a,s,b]=[%.3e,%.3e,%.3e], [ga,gs,gb]=[%.3e,%.3e,%.3e]\n",\
a,s,b,ga,gs,gb);
end
maxerr = 2*tol;
while ( b-a > maxerr ) && nev < maxev,
if gs > ga && gs > gb, # Check assumptions
printf ("brent_line_min : gs > ga && gs > gb\n");
keyboard
end
if ga == gb && gb == gs, break end
# Don't trust poly_2_ex for glued points
# (see test_poly_2_ex).
## if min (b-s, s-a) > 10*seps,
# If s is not glued to a or b and does not
# look linear
## mydet = sum (l([2 3 1]).*f([3 1 2])-l([3 1 2]).*f([2 3 1]))
mydet = sum ([s b a].*[gb ga gs] - [b a s].*[gs gb ga]);
if min (b-s, s-a) > 10*seps && abs (mydet) > 10*seps && \
(t = poly_2_ex ([a,s,b], [ga, gs, gb])) < b && t > a,
# t has already been set
## if t>=b || t<=a,
## printf ("brent_line_min : t is not in ]a,b[\n");
## keyboard
## end
# Otherwise, reduce the biggest of the
# intervals, but not too much
elseif s-a > b-s,
t = max (0.5*(a+s), s-100*seps);
else
t = min (0.5*(s+b), s+100*seps);
end
if abs (t-s) < 0.51*maxerr,
#sayif (verbose, "ungluing t and s\n");
t = s + (1 - 2*(s-a > b-s))*0.49*maxerr ;
end
if a > s || s > b, # Check assumptions
printf ("brent_line_min : a > s || s > b\n");
keyboard
end
xt = args;
args{narg} = x+t*dx;
gt = feval( f, args );
nev++;
if verbose,
printf ("t = %.3e, gt = %.3e\n",t,gt);
end
if t<s, # New point is in ]a,s[
if gt > ga + seps, # Check assumptions
printf ("brent_line_min : gt > ga\n");
keyboard
end
if gt < gs,
b = s; gb = gs;
s = t; gs = gt;
else
a = t; ga = gt;
end
else # New point is in ]s,b[
if gt > gb + seps, # Check assumptions
printf ("brent_line_min : gt > gb\n");
keyboard
end
if gt < gs,
a = s; ga = gs;
s = t; gs = gt;
else
b = t; gb = gt;
end
end
if verbose,
printf ("[a,s,b]=[%.3e,%.3e,%.3e], [ga,gs,gb]=[%.3e,%.3e,%.3e]\n",\
a,s,b,ga,gs,gb);
end
## keyboard
## [b-a, maxerr]
end
s2 = 0.5*(a+b);
args{narg} = x + s2*dx; gs2 = feval (f, args );
nev++;
if gs2 < gs,
s = s2; gs = gs2;
end
if gs > gs0,
printf ("brent_line_min : goes uphill by %8.3e\n",gs-gs0);
keyboard
end## [a, b, ga, gb, nev] = semi_bracket (f, dx, a, narg, args)
##
## Find an interval containing a local minimum of the function
## g : h in reals ---> f (x+h*dx) where x = args{narg}
##
## a < b.
## nev is the number of function evaluations
## Author : Etienne Grossmann <etie...@isr.ist.utl.pt>
## Modified by: Levente Torok <torok...@gmail.com>
## This software is distributed under the terms of the GPL
function [a, b, ga, gb, n] = bracket_min (f, dx, narg, args)
[a,b, ga,gb, n] = semi_bracket (f, dx, 0, narg, args);
if a != 0, return; end
[a2,b2, ga2,gb2, n2] = semi_bracket (f, -dx, 0, narg, args);
n += n2;
if a2 == 0,
a = -b2; ga = gb2;
else
a = -b2; b = -a2; ga = gb2; gb = ga2;
end
## -*- texinfo -*-
## @deftypefn {Function File} {...@var{x0},@var{v},@var{nev}]} cg_min ( @var{f},@var{df},@var{args},@var{ctl} )
## NonLinear Conjugate Gradient method to minimize function @var{f}.
##
## @subheading Arguments
## @itemize @bullet
## @item @var{f} : string : Name of function. Return a real value
## @item @var{df} : string : Name of f's derivative. Returns a (R*C) x 1 vector
## @item @var{args}: cell : Arguments passed to f...@*
## @item @var{ctl} : 5-vec : (Optional) Control variables, described below
## @end itemize
##
## @subheading Returned values
## @itemize @bullet
## @item @var{x0} : matrix : Local minimum of f
## @item @var{v} : real : Value of f in x0
## @item @var{nev} : 1 x 2 : Number of evaluations of f and of df
## @end itemize
##
## @subheading Control Variables
## @itemize @bullet
## @item @var{ctl}(1) : 1 or 2 : Select stopping criterion amongst :
## @item @var{ctl}(1)==0 : Default value
## @item @var{ctl}(1)==1 : Stopping criterion : Stop search when value doesn't
## improve, as tested by @math{ ctl(2) > Deltaf/max(|f(x)|,1) }
## where Deltaf is the decrease in f observed in the last iteration
## (each iteration consists R*C line searches).
## @item @var{ctl}(1)==2 : Stopping criterion : Stop search when updates are small,
## as tested by @math{ ctl(2) > max { dx(i)/max(|x(i)|,1) | i in 1..N }}
## where dx is the change in the x that occured in the last iteration.
## @item @var{ctl}(2) : Threshold used in stopping tests. Default=10*eps
## @item @var{ctl}(2)==0 : Default value
## @item @var{ctl}(3) : Position of the minimized argument in args Default=1
## @item @var{ctl}(3)==0 : Default value
## @item @var{ctl}(4) : Maximum number of function evaluations Default=inf
## @item @var{ctl}(4)==0 : Default value
## @item @var{ctl}(5) : Type of optimization:
## @item @var{ctl}(5)==1 : "Fletcher-Reves" method
## @item @var{ctl}(5)==2 : "Polak-Ribiere" (Default)
## @item @var{ctl}(5)==3 : "Hestenes-Stiefel" method
## @end itemize
##
## @var{ctl} may have length smaller than 4. Default values will be used if ctl is
## not passed or if nan values are given.
## @subheading Example:
##
## function r=df( l ) b=[1;0;-1]; r = -( 2...@{1@} - 2*b + rand(size(l{1}))); endfunction @*
## function r=ff( l ) b=[1;0;-1]; r = (l...@{1@}-b)' * (l...@{1@}-b); endfunction @*
## ll = @{ [10; 2; 3] @}; @*
## ctl(5) = 3; @*
## [x0,v,nev]=cg_min( "ff", "df", ll, ctl ) @*
##
## Comment: In general, BFGS method seems to be better performin in many cases but requires more computation per iteration
## @seealso{ bfgsmin, http://en.wikipedia.org/wiki/Nonlinear_conjugate_gradient }
## @end deftypefn
## Author : Etienne Grossmann <etie...@isr.ist.utl.pt>
## Modified by: Levente Torok <torok...@gmail.com>
## This software is distributed under the terms of the GPL
function [x,v,nev] = cg_min (f, dfn, args, ctl)
verbose = 0;
crit = 1; # Default control variables
tol = 10*eps;
narg = 1;
maxev = inf;
method = 2;
if nargin >= 4, # Read arguments
if !isnan (ctl(1)) && ctl(1) ~= 0, crit = ctl(1); end
if length (ctl)>=2 && !isnan (ctl(2)) && ctl(2) ~= 0, tol = ctl(2); end
if length (ctl)>=3 && !isnan (ctl(3)) && ctl(3) ~= 0, narg = ctl(3); end
if length (ctl)>=4 && !isnan (ctl(4)) && ctl(4) ~= 0, maxev = ctl(4); end
if length (ctl)>=5 && !isnan (ctl(5)) && ctl(5) ~= 0, method= ctl(5); end
end
if iscell (args), # List of arguments
x = args{narg};
else # Single argument
x = args;
args = {args};
end
if narg > length (args), # Check
error ("cg_min : narg==%i, length (args)==%i\n",
narg, length (args));
end
[R, C] = size(x);
N = R*C;
x = reshape (x,N,1) ;
nev = [0, 0];
v = feval (f, args);
nev(1)++;
dxn = lxn = dxn_1 = -feval( dfn, args );
nev(2)++;
done = 0;
## TEMP
## tb = ts = zeros (1,100);
# Control params for line search
ctlb = [10*sqrt(eps), narg, maxev];
if crit == 2, ctlb(1) = tol; end
x0 = x;
v0 = v;
nline = 0;
while nev(1) <= maxev ,
## xprev = x ;
ctlb(3) = maxev - nev(1); # Update # of evals
## wiki alg 4.
[alpha, vnew, nev0] = brent_line_min (f, dxn, args, ctlb);
nev += nev0;
## wiki alg 5.
x = x + alpha * dxn;
if nline >= N,
if crit == 1,
done = tol > (v0 - vnew) / max (1, abs (v0));
else
done = tol > norm ((x-x0)(:));
end
nline = 1;
x0 = x;
v0 = vnew;
else
nline++;
end
if done || nev(1) >= maxev, return end
if vnew > v + eps ,
printf("cg_min: step increased cost function\n");
keyboard
end
# if abs(1-(x-xprev)'*dxn/norm(dxn)/norm(x-xprev))>1000*eps,
# printf("cg_min: step is not in the right direction\n");
# keyboard
# end
# update x at the narg'th position of args cellarray
args{narg} = reshape (x, R, C);
v = feval (f, args);
nev(1)++;
if verbose, printf("cg_min : nev=%4i, v=%8.3g\n",nev(1),v) ; end
## wiki alg 1:
dxn = -feval (dfn, args);
nev(2)++;
# wiki alg 2:
switch method
case 1 # Fletcher-Reenves method
nu = dxn' * dxn;
de = dxn_1' * dxn_1;
case 2 # Polak-Ribiere method
nu = (dxn-dxn_1)' * dxn;
de = dxn_1' * dxn_1;
case 3 # Hestenes-Stiefel method
nu = (dxn-dxn_1)' * dxn;
de = (dxn-dxn_1)' * lxn;
otherwise
error("No method like this");
endswitch
if nu == 0,
return
endif
if de == 0,
error("Numerical instability!");
endif
beta = nu / de;
beta = max( 0, beta );
## wiki alg 3. update dxn, lxn, point
dxn_1 = dxn;
dxn = lxn = dxn_1 + beta*lxn ;
end
if verbose, printf ("cg_min: Too many evaluatiosn!\n"); end
endfunction
## [a, b, ga, gb, nev] = semi_bracket (f, dx, a, narg, args)
##
## Find an interval containing a local minimum of the function
## g : h in [a, inf[ ---> f (x+h*dx) where x = nth (args, narg)
##
## The local minimum may be in a.
## a < b.
## nev is the number of function evaluations.
## Author : Etienne Grossmann <etie...@isr.ist.utl.pt>
## Modified by: Levente Torok <torok...@gmail.com>
## This software is distributed under the terms of the GPL
function [a,b,ga,gb,n] = semi_bracket (f, dx, a, narg, args)
step = 1;
x = nth (args, narg);
args{narg} = x+a*dx; ga = feval (f, args );
b = a + step;
args{narg} = x+b*dx; gb = feval (f, args );
n = 2;
if gb >= ga, return ; end
while 1,
c = b + step;
args{narg} = x+c*dx; gc = feval( f, args );
n++;
if gc >= gb, # ga >= gb <= gc
gb = gc; b = c;
return;
end
step *= 2;
a = b; b = c; ga = gb; gb = gc;
end
# -*- texinfo -*-
#
# @deftypefn {Function File} {...@var{cl}} cl_multinom( @var{x},@var{N},@var{b},@var{calculation_type} ) - Confidence level of multinomial portions
# Returns confidence level of multinomial parameters estimated @math{ p = x / sum x } with confidence interval @var{b}.
# Finite population is also considered.
#
# @subheading Arguments
# @itemize @bullet
# @item @var{x} : int vector : sample frequencies bins
# @item @var{N} : int : Population size that was sampled by x. If N<sum(x), infinite number assumed
# @item @var{b} : real, vector : confidence interval
# if vector, it should be the size of x containing confence interval for each cells
# if scalar, each cell will have the same value of b unless it is zero or -1
# if value is 0, b=.02 is assumed which is standard choice at elections
# otherwise it is calculated in a way that one sample in a cell alteration defines the confidence interval
# @item @var{calculation_type} : string : (Optional), described below
# "bromaghin" (default) - do not change it unless you have a good reason to do so
# "cochran"
# "agresti_cull" this is not exactly the solution at reference given below but an adjustment of the solutions above
# @end itemize
#
# @subheading Returns
# Confidence level.
#
# @subheading Example
# CL = cl_multinom( [27;43;19;11], 10000, 0.05 )
# returns 0.69 confidence level.
#
# @subheading References
#
# "bromaghin" calculation type (default) is based on
# is based on the article
# Jeffrey F. Bromaghin, "Sample Size Determination for Interval Estimation of Multinomial Probabilities", The American Statistician vol 47, 1993, pp 203-206.
#
# "cochran" calculation type
# is based on article
# Robert T. Tortora, "A Note on Sample Size Estimation for Multinomial Populations", The American Statistician, , Vol 32. 1978, pp 100-102.
#
# "agresti_cull" calculation type
# is based on article in which Quesenberry Hurst and Goodman result is combined
# A. Agresti and B.A. Coull, "Approximate is better than \"exact\" for interval estimation of binomial portions", The American Statistician, Vol. 52, 1998, pp 119-126
#
# @end deftypefn
# Copyright (C) 2009 Levente Torok / torok...@gmail.com
# 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 CL = cl_multinom( x, N, b, calculation_type )
k = rows(x);
nn = sum(x);
p = x / nn;
if (nargin < 3)
b = .05;
endif
if (isscalar( b ))
if (b==0) b=0.02; endif
b = ones( rows(x), 1 ) * b;
if (b<0) b=1 ./ max( x, 1 ); endif
endif
bb = b .* b;
if (N==nn)
CL = 1;
return;
endif
if (N<nn)
fpc = 1;
else
fpc = (N-1) / (N-nn); # finite population correction tag
endif
beta = p.*(1-p);
if ( nargin < 4 )
calculation_type = "bromaghin";
endif
switch calculation_type
case {"cochran"}
t = sqrt( fpc * nn * bb ./ beta )
alpha = ( 1 - normcdf( t )) * 2
case {"bromaghin"}
t = sqrt( fpc * (nn * 2 * bb )./ ( beta - 2 * bb + sqrt( beta .* beta - bb .* ( 4*beta - 1 ))) );
alpha = ( 1 - normcdf( t )) * 2;
case {"agresti_cull"}
ts = fpc * nn * bb ./ beta ;
if ( k<=2 )
alpha = 1 - chi2cdf( ts, k-1 ); % adjusted Wilson interval
else
alpha = 1 - chi2cdf( ts/k, 1 ); % Goodman interval with Bonferroni argument
endif
endswitch
CL = 1 - max( alpha );
endfunction
------------------------------------------------------------------------------
Apps built with the Adobe(R) Flex(R) framework and Flex Builder(TM) are
powering Web 2.0 with engaging, cross-platform capabilities. Quickly and
easily build your RIAs with Flex Builder, the Eclipse(TM)based development
software that enables intelligent coding and step-through debugging.
Download the free 60 day trial. http://p.sf.net/sfu/www-adobe-com
_______________________________________________
Octave-dev mailing list
Octave-dev@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/octave-dev
