Nonlinear minimization; group contains NELDER.82P, NELDSORT.82P, TESTNELD.82P
(example), ORBIT.82P (example), LOGISTIC.82P (example), FUNC.82P.
----begin documentation----
This 82u-file was prepared by Mikael Bonnier, http://www.df.lth.se.orbin.se/~mikaelb/.
******************************
* NELDER.82G for the TI-82 *
******************************
ACKNOWLEDGEMENT
NELDER.82G is a group file containing NELDER.82P, NELDSORT.82P, FUNC.82P, and
the examples TESTNELD.82P, LOGISTIC.82P, and ORBIT.82P. They were written by
John Powers and C. B. Wilson of the TI Graphics Team and are released to the
public domain. You may copy and change these programs. This version for the
TI-82 is based on similar programs originally written for the TI-85 (also
located in the SAMPLES directory of this disk). Due to the more limited
variables of the TI-82, the program is less obvious to read, so a mapping
of variable names to the TI-85 program is provided later in this file.
INTRODUCTION
NELDER.82P and its associated subroutine NELDSORT.82P will minimize a
non-linear (or linear) function with one or more parameters, using the
Nelder-Mead simplex method. Nonlinear optimization is a very valuable
tool, especially for modeling and regression. The Nelder-Mead method is
one of the most robust and therefore very suitable for experimentation and
use on a wide range of problems. TESTNELD.82P, LOGISTIC.82P and ORBIT.82P
are application examples for NELDER.82P.
HOW TO USE NELDER
You should provide a program named FUNC (this specific program name is
required) which will compute a real value from your data and the current
value of the parameters to be optimized and place this real value in
variable F. Then execute NELDER, which will try to minimize the quantity
returned in F. The choice of function to minimize is up to you, but for
example if you want least squares minimization you would compute the sum
of the squares of the residuals in program FUNC. You may store data in
variable names of your choice, but the parameters to be optimized in
your model must be in list L1. Also NELDER uses the following variable
names, which must not be used for other purposes:
REALS LISTS MATRICES
N L1 (for coefficients) [D]
J L5 [E]
K L6
L
T
F
In addition to program FUNC, you may want to have a setup program which
sets up data, stopping tolerance, etc., and then calls NELDER. The example
programs TESTNELD.82P, LOGISTIC.82P and ORBIT.82P are of this form.
NELDER uses the value in variable T as a termination criterion. The
iterations stop when the "size" of the simplex (roughly the "distance"
between vertices) is less than T. At each iteration, NELDER displays
the action taken on that iteration (REFLECT, EXPAND, CONTRACT or SHRINK),
the loop count, a list representing the value of F at each vertex, and
the value compared to T for termination. There is no limit on loop
count; you can break the program if it appears it is not going to converge.
The most recent trial for L1 is averaged in [E](1,...) and the last actions
stored in [E](5,...) for "reflection", [E](3,...) for "expansion",
[E](2,...) for "contraction" and [E](4,...) for "shrink". At the beginning
of NELDER is a constant offset of 0.1 stored in L6(10). This value is used
to perturb your initial guess to form the simplex and initialize the
process. If you know your initial guess is much more accurate than 10%,
you can prevent the program "opening up" the scale too much by reducing
this value. If you do this often, you can remove the assignment statement
from NELDER and put it in a calling program.
The example program TESTNELD.82P illustrates a simple use of NELDER. This
example fits a series of Fibonacci numbers to an equation of the form F(n)=
ae^(bn). In the example program L1(1) is "a" and L1(2) is "b", the list L2
is "n" and list L3 contains the exact Fibonacci sequence for the values of
"n" in list L2. The equation for F uses the sum of squared error as
the minimization function. To make the program FUNC more general, the
model equation is stored in Y0.
If you use L2={1,2,3,4,5,6,7,8,9,10}, L3={1,1,2,3,5,8,13,21,34,55} and
T=1e-5, you should get a least-squares error of .12339, L1(1)=.4480817 and
L1(2)=.4810051. This is a good example of why nonlinear minimization is
required for accurate answers to tasks that are often treated in
"linearized" form due to people not having tools like this NELDER program.
For instance, y=ae^(bx) can be linearized by ln(y)=ln(a) + bx and linear
least-square methods used to obtain ln(a) (hence a) and b. However, these
values of a and b are biased by the fact that the sum of squares is
minimized in ln(y) and not y. We can see this effect quickly in the TI-82
by using exponential regression on the two lists above. Since the
exponential regression is y=ab^x (linearized as ln(y)=ln(a) + ln(b)*x), we
just need to take the natural log of b as given by the TI-82 to have the a
and b that we want. For the data above, we obtain a=.488296 and b=.468948.
Substituting these values for L1(1) and L1(2), we can run FUNC and look at
F to see that the sum of squared error is 4.31 (vice 0.12). If we want an
equation for calculating Fibonacci numbers: FIB=round(L1(1) e^(L1(2)*N),0),
the coefficients we obtained using NELDER will return the correct integer
values over the range of our data (and indeed up to N=14), while the
coefficients from the linearized least squares approach only work up to
N=8. In fact, if you try to fit the first 60 Fibonacci numbers with
linearized least squares, your coefficients will result in a sum of squared
error of 1.59e20. Using these coefficients in the equation above only
results in one additional correct Fibonacci number for N=9. But with
NELDER, we can determine the coefficients L1(1)=.44721359550129 and
L1(2)=.48121182505955, which result in a sum of squared error of 0.4407 and
exact results for N from 1 to 60. (Note that to do this we will need to set
the stopping tolerance T to 1E-11.) The lesson here is that for accurate
data modeling to nonlinear equations, you really have to use a nonlinear
minimization tool. While Nelder-Mead is not the fastest algorithm, it is
very robust, doesn't need any derivative information, and isn't bothered too
much by discontinuities or other messy function behavior.
Note: Just as an aside, it may be noted the ratio of successive
Fibonacci numbers approaches the Golden mean (sqrt(5)-1)/2). Based on this
observation, we realize that the nth Fibonacci number should equal
approximately a*b^n; where a = 1/sqrt(5) and b = 2/(sqrt(5)-1).
In the form we fit above therefore, L1(1)=1/sqrt(5) and
L1(2)=ln(2/(sqrt(5)-1)).
A second example LOGISTIC.82P, fits data on the time evolution of an
algae sample taken in the Adriatic Sea to a logistic differential
equation to model the dynamics of a one-species population in an
environment with limited resources. This example is taken from "Fitting
a Logistic Curve to Data" in The College Mathematics Journal, Vol 24, No. 3,
May 1993, pp 247-253. The model equation is m(t) = K / (1 + e^(-r(t - t0)),
where m(t) is the biomass in square millimeters of surface covered by
biomass in a microscope sample versus t in days. The data are:
t 11 15 18 23 26 31 39 44 54 64 74
m(t) .00476 .0105 .0207 .0619 .337 .74 1.7 2.45 3.5 4.5 5.09
We will assign coefficients: L1(1)=K, L1(2)=r, L1(3)=t0. From a statplot,
we note the curve asymtotically approaches an upper bound for large t, so
we might guess K at 6. We also note the curve has an inflection point near
the center of the t range, so we might guess t0 at 40. If we use these
values in the equation at t=74, we can solve for r of about 0.05. This
gives us an initial guess of {6,.05,40}. With T=.001, we get final
coefficients of {5.095,.1213,45.775} and a sum of squared residuals of
.23817.
A third example, ORBIT.82P, fits data on variation in the period of pulsar
PSR1257+12 over time. This data is approximate and obtained from the graph
shown in the announcement article, "A planetary system around the
millisecond pulsar PSR1257+12", Nature, Vol. 355, 9 January 1992, pp
145-147. Another article covering this is "Pulsars, Planets and Pathos" in
Sky & Telescope, May 1992, pp 493-495. The data is contained within the
program itself, as well as initial conditions for the optimization. The
program fits 7 parameters to an equation, which is the sum of two sinusoids
and a constant. Each sinusoid is described by an amplitude, a frequency
and a phase. The time data is in years (actual year - 1990) and the period
data is in nanoseconds (actual period - 6,218,530 ns). This problem is
pretty complex for this method, which is better for five or less
parameters, therefore the calculation time will be fairly long. However,
it is an example of a realistic scientific problem that students have not
generally been equipped to solve. From a sketch of the data points you can
determine possible values for the amplitudes and periods. For instance,
one period seems to be about 0.2 year, so we would use 2*pi/0.2 = 31.4
radians as a possible starting point.
It is more difficult to make a close initial estimate of the phase angles
and if we just start with arbitrary values, it may take several hundred
iterations to get a result with tol=.001 or less. Starting with the
somewhat better initial estimate in the program {.94,35,-2,.67,23,3.2,3.4},
setting the offset L6(10) in NELDER to 0.01 and tol=.01, you can get the
following result in 49 iterations.
First term Second term
Amplitude .93465 .66863
Frequency 34.7047 23.0662
Phase -2.0036 3.2089
Constant 3.3935
Converting the two frequencies into days:
365.25 days * 2pi / 34.7047 = 66.1 days
365.25 days * 2pi / 23.0662 = 99.5 days
As reported in the Nature article, it is believed that two planets are in
orbit around this pulsar, having periods of 66.6 and 98.2 days. With
assumptions about the mass of the pulsar and this data, calculations for
the mass of these planets can also be made.
If you want to see the equation and data together, one way is to do a
Statplot with Xlist of L2 and Ylist of L3. Use ZoomStat to get a window
that covers the data. Recall a copy of Y0 into another Y= equation and
edit the two L2's into X's and graph this equation.
References for the Nelder-Mead technique are:
1. J. A. Nelder and R. Mead, "A Simplex method for function minimization",
Computer Journal, vol. 7, pp 308-313, 1965.
2. W. H. Press, B. P. Flannery, S. A. Teukolsky, W. T. Vetterling,
"Numerical Recipes, The Art of Scientific Computing", Cambridge University
Press, 1987, pp 305-309.
Appendix: Variable name correspondence between NELDER.82P and NELDER.85P.
TI-85 variable name TI-82 variable name
X L1
fv L5
vertex [D]
vavg [E](1,n) where n equals number of coefficients
vcon [E](2,n)
vexp [E](3,n)
vk [E](4,n)
vref [E](5,n)
vtemp [E](6,n)
xi [E](7,n)
swap [E](8,n)
fk L6(1)
cnt L6(2)
test L6(3)
s L6(4)
fref L6(5)
fexp L6(6)
ftemp L6(7)
fcon L6(8)
how L6(9)
offset L6(10)
n9 N
j9 J
k9 K
FX F
----- L (additional loop index)
----begin ascii----
\START82\
\COMMENT=Program file dated 08/13/96, 17:11
\NAME=FUNC
\FILE=FUNC.82P
:sum (\Y0\-\L3\)\^2\\->\F
\STOP82\
\START82\
\COMMENT=Program file dated 08/13/96, 17:11
\NAME=LOGISTIC
\FILE=LOGISTIC.82P
:"\L1\(1)/(1+e^(\(-)\\L1\(2)(\L2\-\L1\(3))))"\->\\Y0\
:{11,15,18,23,26,31,39,44,54,64,74}\->\\L2\
:{.00476,.0105,.0207,.0619,.337,.74,1.7,2.45,3.5,4.5,5.09}\->\\L3\
:{6,.05,40}\->\\L1\
:.001\->\T
:Fix 3
:prgmNELDER
:Float
:Disp \L1\
\STOP82\
\START82\
\COMMENT=Program file dated 08/13/96, 17:11
\NAME=NELDER
\FILE=NELDER.82P
:10\->\dim \L6\
:.1\->\\L6\(10)
:dim \L1\\->\N
:{N+1,N}\->\dim [D]
:N+1\->\dim \L5\
:{8,N}\->\dim [E]
:For(L,1,N)
:\L1\(L)\->\[E](7,L)
:(1-\L6\(10))[E](7,L)\->\\L1\(L)
:\L1\(L)\->\[D](N+1,L)
:End
:prgmFUNC
:F\->\\L5\(N+1)
:For(J,1,N)
:For(L,1,N)
:[E](7,L)\->\\L1\(L)
:End
:If \L1\(J)\<>\0
:Then:(1+\L6\(10))\L1\(J)\->\\L1\(J)
:Else:\L6\(10)\->\\L1\(J)
:End
:For(L,1,N)
:\L1\(L)\->\[D](J,L)
:End
:prgmFUNC
:F\->\\L5\(J)
:End
:prgmNELDSORT
:ClrHome
:1\->\\L6\(2)
:Disp "INITIAL",\L6\(2),\L5\
:Lbl A
:0\->\\L6\(3)
:For(J,2,N+1)
:0\->\\L6\(4)
:For(L,1,N)
:([D](J,L)-[D](1,L))\^2\+\L6\(4)\->\\L6\(4)
:End
:\sqrt\(\L6\(4))\->\\L6\(4)
:max(\L6\(3),\L6\(4))\->\\L6\(3)
:End
:If \L6\(3)\<=\T:Then
:For(L,1,N)
:[D](1,L)\->\\L1\(L)
:End
:Return
:End
:For(K,1,N)
:0\->\\L6\(4)
:For(J,1,N)
:[D](J,K)+\L6\(4)\->\\L6\(4)
:End
:\L6\(4)/N\->\[E](1,K)
:End
:For(L,1,N)
:2[E](1,L)-[D](N+1,L)\->\[E](5,L)
:[E](5,L)\->\\L1\(L)
:End
:prgmFUNC
:F\->\\L6\(5)
:For(L,1,N)
:[E](5,L)\->\[E](4,L)
:End
::\L6\(5)\->\\L6\(1)
:1\->\\L6\(9)
:If \L6\(5)<\L5\(N):Then
:If \L6\(5)<\L5\(1):Then
:For(L,1,N)
:2[E](5,L)-[E](1,L)\->\[E](3,L)
:[E](3,L)\->\\L1\(L)
:End
:prgmFUNC
:F\->\\L6\(6)
:If \L6\(6)<\L5\(1):Then
:For(L,1,N)
:[E](3,L)\->\[E](4,L)
:End
:\L6\(6)\->\\L6\(1)
:2\->\\L6\(9)
:End
:End
:Else
:For(L,1,N)
:[D](N+1,L)\->\[E](6,L)
:End
:\L5\(N+1)\->\\L6\(7)
:If \L6\(5)<\L6\(7):Then
:For(L,1,N)
:[E](5,L)\->\[E](6,L)
:End
:End
:For(L,1,N)
:.5([E](6,L)+[E](1,L))\->\[E](2,L)
:[E](2,L)\->\\L1\(L)
:End
:prgmFUNC
:F\->\\L6\(8)
:If \L6\(8)<\L5\(N):Then
:For(L,1,N)
:[E](2,L)\->\[E](4,L)
:End
::\L6\(8)\->\\L6\(1)
:3\->\\L6\(9)
:Else
:For(J,2,N)
:For(L,1,N)
:([D](1,L)+[D](J,L))/2\->\[D](J,L)
:[D](J,L)\->\\L1\(L)
:End
:prgmFUNC
:F\->\\L5\(J)
:End
:For(L,1,N)
:([D](1,L)+[D](N+1,L))/2\->\[E](4,L)
:[E](4,L)\->\\L1\(L)
:End
:prgmFUNC
:F\->\\L6\(1)
:4\->\\L6\(9)
:End
:End
:For(L,1,N)
:[E](4,L)\->\[D](N+1,L)
:End
:\L6\(1)\->\\L5\(N+1)
:prgmNELDSORT
:\L6\(2)+1\->\\L6\(2)
:If \L6\(9)=1
:Disp "REFLECT"
:If \L6\(9)=2
:Disp "EXPAND"
:If \L6\(9)=3
:Disp "CONTRACT"
:If \L6\(9)=4
:Disp "SHRINK"
:Disp \L6\(2),\L5\
:Disp \L6\(3)
:Goto A
\STOP82\
\START82\
\COMMENT=Program file dated 08/13/96, 17:11
\NAME=NELDSORT
\FILE=NELDSORT.82P
:For(J,1,N)
:For(K,J+1,N+1)
:If \L5\(J)>\L5\(K):Then
:\L5\(J)\->\L:\L5\(K)\->\\L5\(J):L\->\\L5\(K)
:For(L,1,N)
:[D](J,L)\->\[E](8,L):[D](K,L)\->\[D](J,L):[E](8,L)\->\[D](K,L)
:End
:End
:End
:End
\STOP82\
\START82\
\COMMENT=Program file dated 08/13/96, 17:11
\NAME=ORBIT
\FILE=ORBIT.82P
:"\L1\(1)sin (\L1\(2)*\L2\+\L1\(3))+\L1\(4)sin (\L1\(5)*\L2\+\L1\(6))+\L1\\#\
(7)"\->\\Y0\
:{.538,.545,.572,.61,.642,.715,.797,.85,.957,.98,1.215,1.225,1.295,1.30\#\
2,1.308,1.368,1.385,1.4,1.41,1.555,1.652,1.678,1.682,1.695,1.707,1.725,\#\
1.765,1.8,1.85}\->\\L2\
:{2.45,2.19,2.13,2.85,3.4,2.83,4.01,3.45,3.18,4.11,3.96,3.88,3.39,3.54,\#\
3.62,4.28,3.9,3.44,3.12,4.94,2.24,2.58,2.72,3.02,3.48,3.85,3.9,3.4,3.62\#\
}\->\\L3\
:.01\->\T
:{.94,35,\(-)\2,.67,23,3.2,3.4}\->\\L1\
:Fix 3
:prgmNELDER
:Float
:Disp \L1\
\STOP82\
\START82\
\COMMENT=Program file dated 08/13/96, 17:11
\NAME=TESTNELD
\FILE=TESTNELD.82P
:"\L1\(1)e^(\L1\(2)\L2\)"\->\\Y0\
:{1,2,3,4,5,6,7,8,9,10}\->\\L2\
:{1,1,2,3,5,8,13,21,34,55}\->\\L3\
:1\E\\(-)\5\->\T
:{.447,.481}\->\\L1\
:Fix 3
:prgmNELDER
:Float
:Disp \L1\
\STOP82\
----begin uue----
begin 644 NELDER.82G
M*BI423@R*BH:"@!'&7%=`A$-!$8_"P"N``5,
M3T=)4U1)0ZX`K``J70`0,1&#$#%POQ"P70`0,A$070%Q70`0,Q$1$1$J!%X9
M/P@Q,2LQ-2LQ."LR,RLR-BLS,2LS.2LT-"LU-"LV-"LW-`D$70$_"#HP,#0W
M-BLZ,#$P-2LZ,#(P-RLZ,#8Q.2LZ,S,W*SHW-"LQ.C
M70`_"P`B!05.14Q$15(``"(%(`4Q,`2U704_.C$$7040,3`1/[5=``1./PA.
M<#$K3@D$M5P#/TYP,02U700_"#@K3@D$M5P$/]-,*S$K3A$_70`03!$$7`00
M-RM,$3\0,7%=!1`Q,!$17`00-RM,$01=`!!,$3]=`!!,$01<`Q!.<#$K3!$_
MU#]?1E5.0S]&!%T$$$YP,1$_TTHK,2M.$3_33"LQ*TX1/UP$$#7040,A$K700_WET%$#,1/]=!"P"!``5.14Q$4T]25($`?P#3
M2BLQ*TX1/]-+*TIP,2M.<#$1/\Y=!!!*$6Q=!!!+$3[//UT$$$H1!$P^7000
M2Q$$70002A$^3`1=!!!+$3_33"LQ*TX1/UP#$$HK3!$$7`00."M,$3Y<`Q!+
M*TP1!%P#$$HK3!$^7`00."M,$01<`Q!+*TP1/]0_U#_4/]0_"P"B`05/4D))
M5````*(!H`$J70`0,1'"$%T`$#(1@ET!<%T`$#,1$7!=`!`T$<(070`0-1&"
M70%P70`0-A$1<%T`$#<1*@1>&3\(.C4S."LZ-30U*SHU-S(K.C8Q*SHV-#(K
M.C