From: "Eric J. Zmyslowski"
Subject: Curve fitting Program
Here is a curve fitting program that I have written/modified.
Most of it was written by myself, but most of the design and some of the
code was constructed by HP, including {POINT, DEL, BEST, CUF1, CUF2}.
**** This will only work on the Version 2BB. ****
(There probably isn't enough memory in the 28C anyways....)
Put the points on the stack in 1 of 3 forms :
1: (x,y) | a complex number
1: [x,y] | a vector
1: x | a stack order
2: y |
And then run POINT. It will input the points in the statistics matrice.
Be warned that the program requires all data points to be in the domain of
x > 0 , y > 0 . and that Aval > y . This Aval will be understood further on.
To fit a curve to the points you can use the BEST program, this will find the
best fitting curve for those points, or you can run a particular curve formula
to find the best type of that particular curve to fit the points.
BEST will return the values of a, b, and r^2 which is the coefficient
of determination (the correlation).
There are other functions that one can also use:
DEL: this program deletes the most recently added x,y pair.
XPLOT: this program plots the points provided and the most recently
fitted curve.
CSTATS: provides additional statistical information, including an 'F' value
that will better describe the fit of the curve.
(this program is optional, and may be excluded)
ACHNG: this program is used to change the value of Aval.
CLVAR: Purges all unneccasary variables.
X->Y: this is used to compute a Y value for a corresponding X value.
Y->X: this is used to compute an X value for a corresponding Y value.
& will represent the "SIGMA" key.
-> will represent the right arrow key.
SQRT will represent the square root symbol.
--------------------------------------------------
The following equations may be used for fitting data:
LINC: Straight line; y = a * x + b
EXPC: Exponential; y = a * exp(b * x)
PWRC: Power function; y = a * x ^ b
LOGC: Logarithmic; y = a + b * ln(x)
HYPC: Hyperbolic; y = (a * x) / (b + x)
MAXM: Maxima function; y = a * x * exp(b * x)
MINV: Modified inverse; y = a / (b + x)
SIGEXP: Exponential Sigmoid; y = Aval / ( 1 + a * exp(b * x))
SIGMOID: Sigmoid; y = Aval / ( 1 + a * x ^ b)
EXSAT: Exponential saturation; y = a * (1 - exp(b * x))
Where Aval represents a constant.
-------------------------------------------------------------------------------
POINT
{requires ordered pair of numbers in 1 of 3 forms: (x,y); [x,y]; x,y}
<< DUP TYPE
IF DUP 3 == THEN
DROP ARRAY-> DROP
ELSE
IF 1 == THEN
C->R
END
END
DUP IF Aval < THEN
-> x y
<< x y x LN y LN x y /
y INV y x / LN Aval y
- LN Aval y / 1 - LN >>
{ 9 } ->ARRY &+
ELSE
{ 2 } ->ARRY
"Error... Aval < y"
1 DISP
END >>
DEL
{returns an ordered pair in the form [x,y]}
<< &- { 2 } RDM >>
BEST
{returns 3: a
2: b
1: r^2
}
<< 'LINC' DUP EVAL 'PWRC'
CFU2 'LOGC' CFU2 'EXPC'
CFU2 'HYPC' CFU2 'MAXM'
CFU2 'MINV' CFU2 'SIGEXP'
CFU2 'SIGMOID' CFU2
'EXSAT' CFU2 4 DROPN DUP
EVAL #154364d SYSEVAL 5
ROLL 4 DISP >>
XPLOT
<< '&PAR' RCL 1 2 COL&
MAX& MIN& -1 1 FOR x
DUP 1 GET SWAP 2 GET
1 .25 x * + DUP 4 ROLL
* SWAP 3 ROLL * 0 DUP
DUP2 3 DUPN 9 ->ARRY
&+ 2 STEP
SCL& &- &- DROP2 CLLCD
DRW& '&PAR' STO << X
X->Y >> STEQ DRAW EQ >>
CSTATS
{returns the values: variance for n degrees of freedom,
standard error of the estimates,
and the 'F' value.
in the following form: { (variance,n) standarderror Fvalue}
The greater the 'F' value the better the curve fits the data.
}
<< &DAT '&DAT2' STO N& CL&
&DAT2 1 3 ROLL FOR x
DUP x 1 2 ->LIST GET X->Y
DUP SQ 3 PICK x 2 2 ->LIST
GET DUP SQ 4 PICK 3 PICK -
SQ 5 ->ARRY &+ NEXT
TOT SWAP DROP DUP 5 GET N& 2
- DUP 3 ROLL SWAP / DUP 3
ROLL SWAP R->C SWAP SQRT
3 PICK DUP 2 GET SWAP 3 GET
DUP 2 * N& / 6 PICK 1 GET
* 3 ROLL SWAP - SWAP SQ N&
/ + 4 ROLL 5 GET
N& 2 - / IF DUP 0 == THEN
DROP MINR ->NUM
END / 3 ->LIST DUP '&PAR2'
STO &DAT2 '&DAT' STO
'&DAT2' PURGE >>
ACHNG
{changes the value of Aval, and recalculates new stat. if statistic
matrice already exists. requires a real number on level 1.
}
<< DUP 'Aval' STO IFERR '&DAT'
RCL THEN DROP2 ELSE DROP
IF DUP MAX& 2 GET < THEN
"Error... Aval < Y" 1 DISP
ELSE -> a << 1 N& FOR z
'&DAT(z,2)' EVAL DUP a SWAP
- LN '&DAT(z,8)' a 4 ROLL
/ 1 - LN '&DAT(z,9)' STO
STO NEXT >> END END >>
CLVAR
<< { &PAR2 EQ PPAR A B RSQ &PAR
X->Y Y->X &DAT } PURGE >>
Aval
{store any real number in here provided Aval > 0, since this is used by
only two of the curve equations.
}
255
GEQ
<< { "y= Bx+A"
"y= Ae^(Bx)"
"y= Bln(x)+A"
"y= Ax^B"
"y= Ax/(B+x)"
"y= A(1-e^(Bx))"
"y= Aval/(1+Ax^B)"
"y= Aval/(1+Ae(Bx))"
"y= A/(B+x)"
"y= Axe^(Bx)" }
SWAP GET >>
LREG
{ << EXP >>
<< SWAP INV DUP 3
ROLL * SWAP >>
<< SWAP INV DUP 3
ROLL * SWAP >>
<< EXP >>
<< EXP >>
<< EXP >> }
LINC
<< 1 GEQ 2 1
<< B * A + >>
<< A - B / >>
CFU1 >>
EXPC
<< 2 GEQ 4 1
<< B * EXP A * >>
<< A / LN B / >>
CFU1 >>
PWRC
<< 4 GEQ 4 3
<< B ^ A * >>
<< A / B INV ^ >>
CFU1 >>
LOGC
<< 3 GEQ 2 3
<< LN B * A + >>
<< A - B / EXP >>
CFU1 >>
HYPC
<< 5 GEQ 5 1
<< DUP A * SWAP B + / >>
<< DUP B * A 3 ROLL - / >>
CFU1 >>
MAXM
<< 10 GEQ 7 1
<< DUP A * SWAP B * EXP * >>
<< 'Y' STO 'X*EXP(X)=B*Y/A'
'X' { MINR MAXR } ROOT B
/ 'Y' PURGE >>
CFU1 >>
MINV
<< 9 GEQ 6 1
<< B + A SWAP / >>
<< A SWAP / B - >>
CFU1 >>
SIGEXP
<< 8 GEQ 9 1
<< B * EXP A * 1 +
Aval SWAP / >>
<< Aval SWAP / 1 -
A / LN B / >>
CFU1 >>
SIGMOID
<< 7 GEQ 9 3
<< B ^ A * 1 +
Aval SWAP / >>
<< Aval SWAP / 1
- A / B INV ^ >>
CFU1 >>
EXSAT
<< 6 GEQ 8 1
<< B * EXP NEG 1 + A * >>
<< A / NEG 1 + LN B / >>
CFU1 >>
CFU1
<< 'Y->X' STO 'X->Y' STO
OVER COL& LR SWAP ROT
IF DUP 4 >= THEN
3 - LREG SWAP GET
ELSE DROP END
IFERR EVAL THEN DROP2
0 0 0 ELSE SWAP CORR SQ
END 3 DUPN 'RSQ' STO 'B'
STO 'A' STO >>
CFU2
<< DUP EVAL DUP 7 PICK
IF > THEN 5 ->LIST 6
ROLLD 5 DROPN
LIST-> DROP ELSE
5 DROPN END >>
I hope this is of some use to some of you. I know that it has been to me.
This is a powerful program and can easily be modified to include other
curve types.
Eric J. Zmyslowski
Grik@mtus5.bitnet