NELSON%VWSCYG@vmsc.oac.uci.edu (Matthew A. Nelson) (12/19/90)
Below you will find a routine for doing fits to the line Y=a+bX, and
a routine to predict a value of y0 given an x0. I realize that these are
already built into the 48, but the routines below not only provide the
coefficients a, b, and the predicted y0, but also the standard deviations
in these values, Sa, Sb and Sy0. Knowing these errors is often times
crucial, and I have long been amazed that linear regression routines
usually omit them.
The formulae used for the calculations below were gotten from
"Introduction to the Theory of Error", by Yardley Beers (Addison-
Wesley, 1962).
The first routine takes your sigma-dat, and calculates a, Sa, b, Sb
and an estimate of the errors in the y-data, Sy (this is used in the
second routine). It changes the display mode and displays the relevant
information, and also creates the global variables a, Sa, b, Sb and Sy.
%%HP: T(3)A(R)F(.);
\<< RCLF \-> f
\<< 4 FIX \GSY SQ
\GSX^2 * \GSX*Y \GSX \GSY *
* 2 * - N\GS \GSX*Y SQ
* + N\GS \GSX^2 * \GSX SQ
- / \GSY^2 SWAP - N\GS
2 - / \v/ DUP 'Sy'
STO \-> sy
\<< CLLCD N\GS N\GS
\GSX^2 * \GSX SQ - / \v/
sy * DUP 'Sb' STO
\->STR "Sb: " SWAP +
6 DISP N\GS \GSX*Y * \GSX
\GSY * - N\GS \GSX^2 * \GSX
SQ - / DUP 'b' STO
\->STR " b: " SWAP +
5 DISP \GSX^2 N\GS \GSX^2
* \GSX SQ - / \v/ sy *
DUP 'Sa' STO \->STR
"Sa: " SWAP + 4
DISP \GSX^2 \GSY * \GSX
\GSX*Y * - N\GS \GSX^2 *
\GSX SQ - / DUP 'a'
STO \->STR " a: "
SWAP + 3 DISP
"fit to: Y=a+bX" 1
DISP 7 FREEZE
\>> f STOF
\>>
\>>
The second routine takes x0 from the stack and returns the predicted
y0 in level 2, and it's error, Sy0, in level 1. The first routine must
be executed before this one, in order to calculate the required coefficients.
%%HP: T(3)A(R)F(.);
\<< \-> x
\<< x b * a + \GSX^2
\GSX 2 * x * - N\GS x
SQ * + N\GS \GSX^2 * \GSX
SQ - / \v/ Sy *
\>>
\>>
Any comments may be thrown to Matt Nelson (nelson@psroot.ps.uci.edu)