Quadratic Regression and Fitting Data to a Line
emc...@my-deja.com
math classes and beyond. Well, in class today, (I'm the only one with
an HP calculator, because they are the best and other people have never
heard of them) we were discussing quadratic regression and fitting data
to a line. On the HP49G, you can fit statistics to a line, but none of
the fit types have the formula ax to the second power + bx + c, which
is the standard form of quadratic expressions. Does anyone know how to
use these types of equations? I don't have a computer link cable, so I
can't download any programs. Thanks, everyone, for your help.
Sent via Deja.com http://www.deja.com/
Before you buy.
Gilberto E. Urroz
> ....... On the HP49G, you can fit statistics to a line, but none of
> the fit types have the formula ax to the second power + bx + c, which
> is the standard form of quadratic expressions. Does anyone know how to
> use these types of equations? I don't have a computer link cable, so I
> can't download any programs. ...
>
Polynomial Regression==================== Theory
========================================
Note: symbols between parentheses attached to a variable represent
subindices, i.e, a(n) = a, sub-n. Also, the symbol ^ is used to
represent powers, i.e., x-squared is x^2.
Polynomial regression consists of finding the coefficients
b(0), b(1), ...,b(m), that will fit an m-th order polynomial of the form:
y = b(0) + b(1)·x + b(2)·x^2 +... + b(m-1)·x^(m-1) + b(m)·x^m + e
[1]
Here, e represents a random error, which is not typically known.
Typically you will use n data points given as a table:
================
x(1) x(2) ... x(n)
y(1) y(2) ... y(n)
================
If we write the vectors y and b as:
_ _ _ _
| y(1) | | b(0) |:
| y(2) | | b(1) |
| . | and | . |, respectively,
| . | | . |
|_ y(n) _| |_ b(m) _|
and a matrix A given by:
_ _
| 1 x(1) x(1)^2 . . . x(1)^m |
| 1 x(2) x(2)^2 . . . x(2)^m |
| . . . . |
| . . . . |
|_ 1 x(n) x(n)^2 . . . x(n)^m_|,
you can write the multiple-linear regression equation as:
y = A·b + e, [2]
where e is the random error vector,
_ _
| e(1) |
| e(2) |
e = | . | .
| . |
|_ e(n)_|
Because the actual errors are not known, we can write only an
approximation to the actual relationship as (the ^ here is used to
distinguish the original data vector, y, from the approximation
based on the polynomial in [1], y^. I.e., in [3] ^ does not mean a power):
y^ = A·b. [3]
The error vector is therefore:
e = y - y^ = y-A·b [4]
Or, in terms of the vector components:
e(i) = y(i) - y^(i), i = 1,2,...,n. [5]
The sum of squared errors is:
SSE= Sum(i=1 to n of e(i)^2) = eT·e = (y-A·b)T·(y-A·b) [6]
(eT = transpose of e)
The least-square method requires us to find the (m+1) values of b(i),
i = 0,2,..,m, that minimize the value of SSE in [6]. This minimization is
obtained by finding the values of b so that:
d(SSE)/db = 0. [7]
The procedure to calculate the derivatives in [7] are complicated,
however, [7] produces the following matrix equation:
AT·A·b = AT·y, [8] (AT = transpose of A)
which can be solved for b as:
b = (AT·A)^(-1)·AT·y [9]
(For matrices, ^(-1) means the inverse matrix).
This requires you to load the matrix A above with
a(i,1) = 1.0, for i = 1,2,...,n
and
a(i,j) = z(i)^(j-1), for j = 2,3,...,m+1; i = 1,2,...,n.
As before we assume that we have n data points (x(i),y(i)),
i = 1,2,...,n, available in the form of a table. The order of
the polynomial, m, must be provided.
Recall that we cannot, in general, obtain the exact fitting of
equation [10], but an approximation given by:
y^ = b(0) + b(1)·x+ b(2)·x^2 +... + b(m-1)·x^(m-1) + b(m)·z^m
========HOW TO PROGRAM THIS?=================
So, how do you do this in the HP 49 G?
I have developed this program that will do the trick for m = 2,
however, it's an RPN program, so bear with me
(I assume that you are using the algebraic mode):
1. Change to RPN mode: [MODE][+/-][OK]
2. Change the menus from choose boxes to softkey menus
by setting system flag 117, as follows:
a. Check if it is set: [1][1][7][+/-][ENTER][ENTER]
(Make two copies of -117).
Then, enter: [CAT][ALPHA][ALPHA][F][S][ALPHA]
(This will select the command FS? from the catalog)
Press [OK].
If you get a 1, the flag is set and you need not do anything else.
Go to step 3, below.
b. If you get a 0, then set the flag by using:
[CAT][ALPHA][ALPHA][S][F][ALPHA][OK]
or just type [ALPHA][ALPHA][S][F][ENTER]
without using the catalog.
3. Type in the following program: (-> stands for the arrow key,
i.e., [rs][0]. Some commands that use the arrow as part of
their names are: OBJ->, ->ARRY, COL->. The left and right-
shift keys are represented by: ls = left shift[blue],
rs = right shift[RED]).
Some notes while typing the program:
* The program symbol << >> is generated by pressing [rs][ + ],
it provides the closing quotes, so you only have to type it at the
beginning of each set of << in the listing below. Same with
[rs][ ' ],
which generates the single quotes ' '. To get out of the quotes
you
will have to use the right arrow key ( from the upper right
corner
of the keyboard.
* The arrow -> is obtained by using: [rs][0],
* To type lower case letters: [ALPHA][ls][(letter)],
e.g., to type j: [ALPHA][ls][J]
* SQ is the [ls][x^2] (x squared function)
* To get the commands OBJ-> and ->ARRY use:
[ls][PRG][TYPE], the programs will be available as softkey menus.
* To get the command SIZE: [ls][PRG][LIST][ELEM][SIZE]
* To get FOR and NEXT: [ls][PRG][BRCH][ELEM][FOR] the programs will
be
available as softkey menus.
* To get TRAN, use: [ls][MATRICES][OPER][NXT][NXT][TRAN].
* To get COL->, use: [ls][MATRICES][CREAT][COL][COL->]
* To get INV, use: [1/x]
* To get DUP, SWAP, DROP, use: [ls][PRG][STACK], the programs will be
available as softkey menus.
* Here is the program listing:
<< -> xx yy
<< xx SQ -> x2
<< xx SIZE 1 SWAP
FOR j 1
NEXT xx SIZE
->ARRY xx OBJ-> ->ARRY 3 COL->
DUP TRAN -> A AT
<< AT A * INV AT
* yy OBJ-> ->ARRY *
>> OBJ-> DROP 'x^2'
* SWAP 'x' * +
SWAP + 'y' SWAP =
>>
>>
>>
The last thing you type is the = sign above, since the >> closing quotes
for
programs are already present. When done, press [ENTER].
4. Save the program in an appropriate variable, say, FITSQ. (Note: to get
back your list of variables after typing, press [VAR]). Use:
[rs]['][ALPHA][ALPHA][F][I][T][S][Q][ENTER] [STO]
5. The program is ready.
=========HOW TO USE THE PROGRAM?=================
Suppose you have a table that you want to fit to
a quadratic equation:
===============================
x 0 0.5 1.0 1.5 2.0 2.5
y 2.8 1.1 0.1 0.2 1.2 2.9
===============================
Thus, n = 6, and, of course, m = 2.
WHILE STILL IN RPN MODE:
Enter a list of data for x and one for y, then press the softkey for FITSQ.
For the problem above enter the following
(Note that elements of a list in RPN mode are separated by spaces):
[ls][{}][0][spc][.][5][spc][1][spc][1][.5]
[spc][2][spc][2][.][5][ENTER]
[ls][{}][2][.][8][spc][1][.][1][spc][.][1]
[spc][.][2][spc][1][.][2][spc][2][.][9]
[ENTER]
[FISTQ]
The result will be:
'y=1.7857*x^2.+-4.4128*x+2.8071'
(Assuming 4 decimals in the output)
=========HOW TO USE THE PROGRAM IN ALGEBRAIC MODE?
Change back to algebraic mode (assuming you were in that awful mode from
the
beginning): [MODE][+/-][OK]
I assume that the program [FITSQ] is visible in the menu, otherwise press
[VAR] to
recover your variable list.
Press the softkey corresponding to FITSQ, then (without pressing [ENTER])
type
[ls][( )]
[ls][{}][0][rs][,][.][5][rs][,][1][rs][,]
[1][.5][rs][,][2][rs][,][2][.][5][right arrow key]
[rs][ , ]
[ls][{}][2][.][8][rs][,][1][.][1][rs][,][.][1][rs][,]
[.][2][rs][,][1][.][2][rs][,][2][.][9][ENTER]
(Note: elements of a list in ALGEBRAIC mode must be separated by commas)
You should get the same result.
=========================THE END======================
I hope this helps. And, I agree with you: the HP is the best.
Greetings.
Gilberto E. Urroz
Department of Civil and Environmental Engineering
Utah State University
Logan,Utah
Keith J. Farmer
news:878d1k$mk1$1...@nnrp1.deja.com...
> Below is the polynomial regression program PREGR
> [Copyright (c) 1992 Hewlett-Packard Company]
Yeah, the one I made just calculates the coefficients, and doesn't draw the
plot. It returns, actually, the row-reduced matrix solution to avoid
problems with ill-conditioned situations.
I suppose I should post the update to that library, fixing the problem with
singular matrices, no?
John H Meyers
> I'm the only one [in class] with an HP calculator, because
> they are the best and other people have never heard of them
That's what's wrong with the educational system :)
Try to acquaint the teachers with it; they need an upgrade ;)
> we were discussing quadratic regression
> On the HP49G, you can fit statistics to a line, but
> none of the fit types have the formula ax^2 + bx + c
Do all the other calculators have quadratic regression built in?
> I don't have a computer link cable, so I can't download any programs.
If you use your calculator a lot,
it may be worthwhile to get or make a cable, e.g. see
<http://www.muffet.com/samsoncables/hp/pc49cable.html>
or <http://hpcalc.org/buying.html> for a shopping guide.
-------------------------------------------------------
Below is the polynomial regression program PREGR
[Copyright (c) 1992 Hewlett-Packard Company]
from <http://www.hpcalc.org/math/numeric/pregr.zip>
or <ftp://ftp.hp.com/pub/calculators/hp48g/math/numeric/pregr.zip>
(because it is hard to read from that file);
everything below is quoted from that source.
PREGR takes [x, y] data from the two-dimensional reserved variable
{Sigma}DAT and returns a polynomial least square fit to that data.
The order of the polynomial is determined by a positive, non-zero,
integer in level 1 of the stack. Given an integer of 4,
the program would return a polynomial of the form:
'Y = a*X^4 + b*X^3 + c*X^2 + d*X + e'
where a, b, c, d, and e would be real coefficients
of the independent variable X.
The program first plots a scatter diagram of the data.
It then determines the coefficients and defines the equation.
It then plots the equation against the scatter diagram.
After plotting the equation the routine terminates in the SOLVER
with the equation in level 1 of the stack and in
the reserved variable EQ. You may use the solver
to compute Y given X or X given Y.
No error checking of input data is performed. The x data need not be
evenly spaced as is required in some polynomial regression programs.
This program is based on an HP-28S program submitted to HP
by Dr. Robert C. Wyckoff.
Example
-------
Fit the data below to a polynomial of degree three.
x: 1.00 1.25 1.50 1.75 2.00 2.25 2.50 2.75 3.00
y: 2.72 3.49 4.48 5.75 7.39 9.49 12.1 15.6 20.00
Keystrokes: Comments:
[Set Approximate mode first, on the HP49]
97 MENU Bring up statistics menu.
(type the above and press ENTER)
|CL{Sigma}| Clear {Sigma}DAT.
[1 2.72 Key in first [x,y] pair as vector.
|{SIGMA}+| Put first pair in {Sigma}DAT.
1.25 [ENTER] 3.49 |{SIGMA}+| Input next data point.
.
.
. Key in rest of the data.
.
3 [ENTER] 20 |{SIGMA}+| Key in last data point.
[VAR] Select variable menu containing |PREGR|.
3 Define the degree as 3.
|PREGR| See a scatter diagram.
See a plot of the fit.
See the equation in the SOLVE application
In the program text below, \GS represents the capital Greek Sigma
(summation sign, red S on HP49), \-> represents the
right-arrow character (red zero on HP49),
and \<< \>> represent program delimiters (red + on HP49);
if you enter the program after setting Approximate mode
on the HP49, you will not need to type any decimal points.
\<< N\GS \-> ORD n \<< SCATRPLOT 0. ORD FOR J
'\GS(I=1.,n,\GSDAT(I,2.)*\GSDAT(I,1.)^J)' EVAL
NEXT { 'ORD+1.' } \->ARRY N\GS 1. ORD
FOR J '\GS(I=1.,n,\GSDAT(I,1.)^J)' EVAL
NEXT ORD 1. + ORD 2. *
FOR J ORD DUPN '\GS(I=1.,n,\GSDAT(I,1.)^J)' EVAL
NEXT { 'ORD+1.' 'ORD+1.' } \->ARRY / -3. CF OBJ\-> OBJ\->
DROP 1. - PATH DUP OBJ\-> 1. SWAP
START EVAL 'X' PURGE NEXT EVAL \-> ORD
\<< 0. SWAP 0. ORD 1. - FOR J 'X' ORD J
- ^ * + SWAP NEXT + 'Y' SWAP = \>>
\>> DUP STEQ 'X' INDEP FUNCTION DRAW TEXT 30. MENU \>>
-----------------------------------------------------------
With best wishes from: John H Meyers <jhme...@mum.edu>
Keith Farmer
answered...
IIRC, the least-squares fit for a jth-degree polynomial to a set of data is
[ Sum(x^0) Sum(x^1) Sum(x^2) ... Sum(x^j) ] [ a0 ] [ Sum(y) ]
[ Sum(x^1) Sum(x^2) Sum(x^3) ... Sum(x^(j+1)) ] [ a1 ] [ Sum(y*x^1) ]
[ Sum(x^2) ... ... ... ... ] * [ a2 ] = [ Sum(y*x^2) ]
[ ... ... ... ... ... ] [ .. ] [ ... ]
[ Sum(x^j) ... ... ... Sum(x^(2j)) ] [ aj ] [ Sum(y*x^j) ]
In the interests of avoiding round-off error as well as to speed calculation
time, I'd like to transform this into integrals, rather than summations. Is
it valid to simply replace all the Sum with Int, perform the integrals, and
solve the resulting system?
In particular, x=0..1, and the first matrix under integrals would then be
invertable (checked w/ RREF, the elements are a.mn=1/(m+n-1)), and the
matrix on the right side has y=Sum(a.i*cos(b.i+c.i*x)), for which y*x^n
(n=0, 1, 2, ...) has closed-form integrals (checked with Mathematica).
I'm interested in making sure this can happen, before I attempt to download
the rather substantial data set which contains the a,b,c used to create
y(x).