#!/bin/sh
# This is a shell archive (produced by GNU sharutils 4.2).
# To extract the files from this archive, save it to some FILE, remove
# everything before the `!/bin/sh' line above, then type `sh FILE'.
#
# Made on 2001-07-05 17:10 CDT by .
# Source directory was `/people/biostat2/therneau/ibm'.
#
# Existing files will *not* be overwritten unless `-c' is specified.
# This format requires very little intelligence at unshar time.
# "if test", "echo", "mkdir", and "sed" may be needed.
#
# This shar contains:
# length mode name
# ------ ---------- ------------------------------------------
# 2728 -rw-r--r-- Readme
# 4060 -rw-r--r-- simplex.d
# 3575 -rw-r--r-- simplex.s
#
echo=echo
if mkdir _sh23484; then
$echo 'x -' 'creating lock directory'
else
$echo 'failed to create lock directory'
exit 1
fi
# ============= Readme ==============
if test -f 'Readme' && test "$first_param" != -c; then
$echo 'x -' SKIPPING 'Readme' '(file already exists)'
else
$echo 'x -' extracting 'Readme' '(text)'
sed 's/^X//' << 'SHAR_EOF' > 'Readme' &&
X This implements the Nelder-Mead simplex algorithm for function
Xminimization. The code is entirely in S-plus.
X
X The algorithm is based on a very simple geometric idea. If there
Xare d parameters, the program will maintain a list of d+1 points, each
Xa guess at the solution. These points can be thought of as enclosing
Xa triangular volume, or simplex, in d dimensional space. At each
Xiteration the algorithm attempts to replace the worst of its guesses
Xwith something better by doing a reflection step: in 3-space imagine
Xthe tetrahedron of the four current guesses as a sock with the 3
X"better" guesses ringing the opening and the bad guess at the toe,
Xthen turn the sock inside out by pulling the toe through. Other step
Xtypes ajdust the step length by shrinking or expanding the point cloud
Xof guesses. At the end of iteration, the points should enclose the
Xtrue solution; the point with the smallest function value is returned
Xas the answer. The algorithm is invariant to rescaling of the
Xparameters, i.e., replacement of p1 with 100* p1, where p1 is the first
Xparameter.
X
XThe basic "run away from the worst" thrust of the algorithm makes it
Xvery robust against being trapped in local minima, see the "smooth +
Xwobble" function in the help page for an example. The program also avoids
Xunwarrantted huge steps that sometimes cause the more common
X"extrapolate towards the best" methods to fail. A reasonable starting
Xstep size (for the author) has been to make the interval (start,
Xstart+step) for each parameter be such that I would a priori believe
Xthat there is about a 50% chance that the final solution for that
Xparameter is in the interval.
X
XA downside to the algorihm is that once it gets in the vicinity of the
Xfinal solution it can take nearly d iterations to shrink each
Xdimension of the solution region by half, far different from the
Xsuper-linear or quadratic behavior of other methods. Do not be
Xsurprized if reduction of the tolerance greatly increases the number
Xof iterations. The author has had success using
Xthis method to generate good starting values for other techniques.
X
X---
X
X This code may be copied or redistributed without restriction.
X
X Terry Therneau
X therneau.terry@mayo.edu
X
X
X----
X A bug in the code was fixed 7/2001, with the initial values. In
X3 dimensions, for example, if the starting estimates were (a,b,c) with
Xinitial step sizes of (1,2,3), the four starting points for the
Xsimplex should be:
X (a,b,c) (a+1, b, c) (a, b+2, c) and (a, b, c+3)
X
XDue to an inadvertent matrix transpose the routine was using
X (a,b,c) (a+1, a, a) (b, b+2, b) and (c, c, c+3)
X
X The only positive note to this is that the minimizer is robust enough
Xso that all my test examples still worked.
SHAR_EOF
: || $echo 'restore of' 'Readme' 'failed'
fi
# ============= simplex.d ==============
if test -f 'simplex.d' && test "$first_param" != -c; then
$echo 'x -' SKIPPING 'simplex.d' '(file already exists)'
else
$echo 'x -' extracting 'simplex.d' '(text)'
sed 's/^X//' << 'SHAR_EOF' > 'simplex.d' &&
X.BG
X.FN simplex
X.TL
XMinimize a function using the Nelder-Mead simplex algorithm
X.DN
XThis will minimize an arbitrary S function with respect to a set
Xof parameters. No deriviatives are required.
X.CS
Xsimplex(func, start, step, ..., tol=0.001, maxiter=1000, trace=F, verbose=F)
X.RA
X.AG func
Xthe function to be minimized. It must have a vector of parameters
Xas its first argument, and return a single value.
X.AG start
Xvector of starting estimates for the parameters.
X.AG step
Xinitial step size for each parameter.
X.OA
X.AG ...
Xany further arguments to the called function. They
Xwill be passed through unchanged as arguments 2, 3, ... of func.
X.AG tol
Xrelative convergence tolerance.
XIteration stops when the relative range of function values within the
Xsolution's current solution set is < tol, i.e.,
X|best-worst| < tol * mean(|best|, |worst|).
X.AG maxiter
Xmaximum number of iterations
X.AG trace
Xprint out a trace of the iteration path, i.e., the solution points
Xthat are attempted.
X.AG verbose
Xinclude the complete set of attempted solutions and their function values
Xin the output.
X.RT
Xa list with elements
X.AG result
Xthe final, minimized value of the function
X.AG x
Xthe parameter values of the solution
X.AG iter
Xthe number of iterations performed
X.AG stepcount
Xa vector containing the number of iteration steps of each type ---
Xreflection, expansion, contraction and shrinkage --- performed in
Xfinding the solution.
XThis is mostly of acedemic interest.
X.AG trace
Xoptional, a matrix containing a trace of the iteration path.
X.DT
X The Nelder-Mead algorithm is based on a very simple geometric idea.
XIf there
Xare d parameters, the program will maintain a list of d+1 points,
Xeach a guess at the solution.
XThese points can be thought of as enclosing a triangular volume, or
Xsimplex, in d dimensional space. At each iteration the algorithm
Xattempts to replace the worst of its guesses with something better
Xby doing a reflection step: in 3-space imagine the tetrahedron of the
Xfour current guesses as a sock with the 3 "better" guesses ringing
Xthe opening and the bad guess at the toe, then turn the sock inside
Xout by pulling the toe through. Other step types ajdust the step length
Xby shrinking or expanding the point cloud of guesses.
XAt the end of iteration, the points should enclose the true solution;
Xthe point with the smallest function value is returned as the answer.
X.PP
XThe basic "run away from the worst" thrust of the algorithm makes it
Xvery robust against being trapped in local minima,
Xand to the unwarrantted huge steps that sometimes cause the
Xmore common "extrapolate towards
Xthe best" methods to fail.
XA reasonable starting step size (for the author)
Xhas been to make the interval (start, start+step) for each parameter
Xbe such that I would a priori believe that there is about a 50%
Xchance that the final solution for that parameter is in the interval.
X.PP
XA downside to the algorihm is that once it gets in the vicinity of the
Xfinal solution it can take nearly d iterations to shrink each dimension
Xof the solution region by half, far different from the super-linear or
Xquadratic behavior of other methods.
XDo not be surprized if reduction of the tolerance greatly increases the
Xnumber of iterations.
X.SH REFERENCES
XNelder, J.A. and Mead, R. 1965, Computer Journal, vol 7, p 308-313.
X
XPress, W.H., Teukolsky, S.A., Vetterling, W.T. and Flannery, B.P. (1992).
XNumerical Recipies in C, second edition. Cambridge University Press.
X.EX
X# Minimize a bouncy smooth function in 2 dimensions
X# z = log(1 +x^2 + y^2) + sin(8x)/4
X#
X> wobble <- function(coef) {
X x <- coef[1]
X y <- coef[2]
X log(1 +x^2 + y^2) + sin(8*x)/4
X }
X> nfit <- simplex(wobble, start=c(-2,-1), step=c(3,3))
X
X#plot a cross section of the function: most minimizers would have trouble
X> x <- seq(-5,5,length=150)
X> plot(x, log(1+x^2) + sin(8*x)/4)
X
X# SLOW least squares
X> rss <- function(beta, x, y) {
X resid <- y - x%*% beta
X var(resid)
X }
X> lfit <- simplex(rss, rep(-1,5), rep(2,5), tol=1e-5,
X x=swiss.x, y=swiss.fertility)
X.KW minimization
X.WR
SHAR_EOF
: || $echo 'restore of' 'simplex.d' 'failed'
fi
# ============= simplex.s ==============
if test -f 'simplex.s' && test "$first_param" != -c; then
$echo 'x -' SKIPPING 'simplex.s' '(file already exists)'
else
$echo 'x -' extracting 'simplex.s' '(text)'
sed 's/^X//' << 'SHAR_EOF' > 'simplex.s' &&
X#
X# S function to minimize using the simplex method
X# slow but sure for bumpy functions
X#
X# Trace arg: print out results as I go
X# verbose arg: return a list of my successful steps
Xsimplex <- function(func, start, step, ...,
X tol=.001, maxiter=1000, trace=F, verbose=F) {
X #
X # Initialize
X #
X nparm <- length(start)
X if (length(step) != nparm) stop ("Wrong length for step vector")
X points <- outer(start, rep(1,nparm))
X points <- cbind(start, points + diag(step))
X
X result <- double(nparm+1)
X for (i in 1:(nparm+1)) result[i] <- func(points[,i], ...)
X
X ord <- order(result)
X points <- points[,ord]
X result <- result[ord]
X
X step <- 'reflect'
X if (verbose) {
X tpoints <- rbind(points, result)
X if (trace) print(tpoints)
X }
X
X stepcount <- rep(0,4)
X history <- integer(maxiter)
X for (iter in 1:maxiter) {
X ii <- iter
X temp <- range(result)
X if ( diff(temp) < (mean(abs(temp)) *tol)) break
X
X if (step=='reflect') {
X history[iter] <- 1
X # Walk away from the worst (last on list)
X center <- apply(points[,1:nparm], 1,mean) #center of opposite face
X newx <- 2*center - points[,nparm+1]
X newy <- func(newx, ...)
X if (trace) cat("r", format(newx), format(newy), "\n")
X
X if (newy < result[1]) {
X step <- 'expand'
X points <- cbind(newx, points[,1:nparm])
X result <- c(newy, result[1:nparm])
X if (verbose) tpoints <- cbind(tpoints, c(newx, newy))
X }
X else if (newy > result[nparm]) {
X step <- 'contract'
X if (newy > result[nparm+1]) newx <- points[,nparm+1]
X }
X else {
X ord <- sum(newy > result)
X points <- cbind(points[,1:ord], newx, points[,(ord+1):nparm])
X result <- c(result[1:ord], newy, result[(ord+1):nparm])
X if (verbose) tpoints <- cbind(tpoints, c(newx, newy))
X }
X }
X else if (step=='expand') {
X # Last step worked great: go further
X history[iter] <- 2
X newx <- 2*newx - center
X newy <- func(newx, ...)
X if (trace) cat("e", format(newx), format(newy), "\n")
X
X if (newy < result[1]) { #keep going
X points[,1] <- newx
X result[1] <- newy
X if (verbose) tpoints <- cbind(tpoints, c(newx, newy))
X }
X else step <- 'reflect'
X }
X
X else if (step=='contract') {
X # Last attempt got worse, try to recover
X history[iter] <- 3
X newx <- .5*(center + newx)
X newy <- func(newx, ...)
X if (trace) cat("c", format(newx), format(newy), "\n")
X
X if (newy > result[nparm]) {
X # still the worst, do a shrinkage towards best point
X #
X stepcount[4] <- stepcount[4] +1
X points[,-1] <- 0.5*(points[,1] + points[,-1])
X for (i in 1+(1:nparm))
X result[i] <- func(points[,i], ...)
X ord <- order(result)
X points <- points[,ord]
X result <- result[ord]
X if (trace) {
X cat ("s\n")
X print(rbind(points, result))
X }
X }
X else if (newy < result[1]) {
X points <- cbind(newx, points[,1:nparm])
X result <- c(newy, result[1:nparm])
X if (verbose) tpoints <- cbind(tpoints, c(newx, newy))
X }
X else {
X ord <- sum(newy > result)
X points <- cbind(points[,1:ord], newx, points[,(ord+1):nparm])
X result <- c(result[1:ord], newy, result[(ord+1):nparm])
X if (verbose) tpoints <- cbind(tpoints, c(newx, newy))
X }
X step <- 'reflect'
X }
X }
X
X stepcount[1:3] <- c(sum(history==1), sum(history==2), sum(history==3))
X names(stepcount) <- c("reflect", "expand", "contract", "shrink")
X if (verbose) {
X list(result=result[1], x=points[,1], iter=ii, stepcount=stepcount,
X trace=t(tpoints))
X }
X else list(result=result[1], x=points[,1], iter=ii, stepcount=stepcount)
X }
X
X
X
SHAR_EOF
: || $echo 'restore of' 'simplex.s' 'failed'
fi
rm -fr _sh23484
exit 0