fx-5800p: Least Squares Fitting

fx-5800p: Least Squares Fitting

Program LSQ (least squares) for the Casio fx-5800p. The program LSQ (least squares) allows us to fit a polynomial or a multiple-linear regression line given certain data and outcomes.

Variables:

Mat A = independent data matrix
Mat B = dependent data matrix, 1-column matrix
Mat C = coefficient matrix. This is the matrix LSQ solves for.
Mat D = predictive values

Calculations:

Coefficients
Mat C = (Mat A^T * Mat A)^-1 * Mat A^T * Mat B = (Mat A)^+ * Mat B

Predictive Values
Mat D = Mat A * Mat C

Prog "LSQ"
"LSQ"
"A = IND. DATA"? → Mat A
"B = DEP. DATA"? → Mat B
( Trn( Mat A ) * Mat A )^-1 * Trn( Mat A ) * Mat B → Mat C
"C = COEFS."
Mat C ◢
Mat A * Mat C → Mat D
"D = PRED."
Mat D

Examples

Multiple Linear Regression

Fit the following data to y = a * x1 + b * x2 + c:
(x1, x2, y):
(1, .6, .45)
(2, .3, .49)
(3, .2, .36)
(4, .8, .36)
(5, .6, .39)

Note that the equation has a constant term c. Set up our matrix A with three columns: one corresponding to the term a * x1, the second for the term b * x2, and a column of ones for the constant column c.

Mat A =
[[ 1, .6, 1]
[2, .3, 1]
[3, .2, 1]
[4, .8, 1]
[5, .6, 1]]

Matrix B will have the dependent values, in this case a 5 x 1 matrix.

Mat B =
[[ .45 ]
[ .49 ]
[ .36 ]
[ .38 ]
[ .39 ]]

With the matrices set up, run LSQ. You can enter the required matrices manually (with [ ALPHA ] [ ln ] for the left bracket, [ ALPHA ] [ x^ ] for the right bracket) or used pre-stored matrices. The fx-5800p has room for only six matrices (Mat A - Mat F).

Results:
Coefficients, Mat B =
[[ -0.02381395349 ]
[ 0.0162790697 ]
[ 0.4773023256 ]]

This corresponds to the equation for the line:

y = -0.02381395349 * x1 + 0.0162790697 * x2 + 0.4773023256

Predicted y values, Mat C =
[[ 0.463255813951]
[ 0.4345581395 ]
[ 0.4091162791 ]
[ 0.3950697674 ]
[ 0.368 ]]

Polynomial Regression

Fit the quadratic polynomial of form y = a*x^2 + b*x + c with data:
(-3, -2.25)
(-1, 0)
(1, -0.11)
(3, 2.23)
(5, 5.24)

We have one variable. However the equation has both x and x^2. Hence, we will set up matrix A as follows: first column has x^2 values, second has x values, and third column contains ones because it corresponds to the constant term (c).

Mat A =
[[ 9, -3, 1]
[1, -1, 0]
[1, 1, -0.11]
[9, 3, 1]
[25, 5, 1]]

Again, matrix B has the dependent values:

Mat B =
[[ -2.24 ]
[ 0 ]
[ -0.11 ]
[ 2.23 ]
[ 5.24 ]]

Results:

Mat C =
[[ 0.07125 ]
[ 0.717 ]
[ -0.33425 ]]

which corresponds to: y = 0.07125 * x^2 + 0.717 * x - 0.33425

Mat D =
[[ -1.844 ]
[ -0.98 ]
[ 0.454 ]
[ 2.458 ]
[ 5.032 ]]

Even though the least-squares method isn't perfect, it is super powerful. Until next time, have a great day, stay healthy, and stave off any viruses because they have gone crazy this winter!

Eddie

Source: http://en.m.wikipedia.org/wiki/Regression_analysis

This blog is property of Edward Shore. 2015

Length of a Polynomial Segment - Part 1; and the HP 50g Program MPFIT

The next two blog entries are going to deal with finding the length of a polynomial segment.

Polynomial Segments

What is a polynomial segment? I could not find a formal definition, so I'll describe it is this: similar to a line segment where a pair of points are connected with a line, a polynomial segment is a set of points that are connected by polynomial curve. The degree of the polynomial is determined by the number of points, less one. For instance, if we have 3 points, the polynomial that connects them the order 2 (quadratic). A cubic polynomial (order 3) is used to connect a set of 4 points.

The Goal

We have a set of points {(x0, y0), (x1, y1), (x2, y2),... (xn, yn)}. Connect the points using a polynomial. What is the length of that segment?

General Method

1. Find the polynomial to connect the set of points. Order of the elements is important. This is determined by using the matrix calculation (X^T X)^-1 (X^T y) to determine the coefficients. The resulting matrix contains the coefficients of the polynomial; with each term being a coefficient of x of increasing power, from 0 to n-1.

2. Use the arc length formula ∫ ( √ (1 + (f'(x))^2 ) dx, a, b) where a is the minimum of the x values and b is the maximum of the x values.

We recommend using a calculator or math software that can handle matrices and integrals for this procedure.

But what is X and y?

y represents a vector of y-values of the set of points.

X is a Vandermonde Matrix. A Vandermonde matrix is a n × n matrix is formed with each row containing a geometric progression of each term. The first column, each term is raised to the 0th power, the second column each term to the 1st, the third column each term is raised to the second power, and so on until the nth row where each term is raised to the n-1 power.

For instance, a Vandermonde matrix consisting of the set [a, b, c, d] would look like this:

Example 1

Find a polynomial segment that fits the points (0,0), (4,8), and (8,6). Find the length of that segment. (See the diagram above).

The x values are [0, 4, 8] with xmin = 0 and xmax = 8. The y values are [0, 8, 6]. Let X be the Vandermonde matrix with the set [0, 4, 8] and the vector y be represented by the corresponding y-values.

Now we need to calculate (X^T X)^-1 (X^T y).

Then:

The resulting matrix contains the coefficients of the required polynomial. Since the vector has 3 entries, the order of the polynomial is 3-1=2.

Going top-down:
0 is the constant (coefficient of x^0),
3.25 is the coefficient of x,
and -0.3125 is the coefficient of x^2.

Our polynomial is f(x) = 3.25x - 0.03125x^2.

To find the length, use the arc length formula using end points xmin = 0 and xmax = 8. We find that the length of the polynomial segment connecting (0,0), (4,8), and (8,6) is about 14.23920.

Example 2

What is the length of the polynomial segment connected with the points (-2,5), (0,-6), (3,2), and (4,1). We are working with 4 points, the required polynomial has the order 4-1=3. (Cubic polynomial)

The Vandermonde matrix X with the vector y are as follows:

Calculating (X^T X)^-1 (X^T y) yields (approximately):

[ [ -6 ]
[ 2.98333 ]
[ 2.725 ]
[ -0.75833 ] ]

To form the cubic polynomial:
-6 is the constant,
2.98333 is the coefficient of x,
2.725 is the coefficient of x^2, and
-0.75833 is the coefficient of x^3.

With xmin = -2 and xmax = 4, we find the length of this segment is approximately 32.60952. (See below).

HP50g Program MPFIT

The HP 50g program MPFIT fits a polynomial to the set of points {(x0, y0), (x1, y1), (x2, y2), ... , (x_(n-1), y_(n-1))}.

Input:
2: vector of x coordinates
1: vector of y coordinates

Output:
2: Coefficients in matrix form [ [a0] [a1] [a2] ... [a_(n-1)] ]
1: Polynomial of X in the form a0 + a1 * X + a2 * X^2 + ... + a_(n-1) * X^(n-1)

Program MPFIT (size 183.5 bytes)

<< DUP SIZE OBJ→ DROP → MX MY N
<< 'X' PURGE
MY OBJ→ 1 + →ARRY
MX VANDERMONDE
LSQ → MR
<< MR DUP
{1, 1} GET
2 N FOR K
MR K 1 2 →LIST GET
'X' K 1 - ^ * +
NEXT >> >> >>

Example:
Input:
2: [1, 3, 6]
1: [6, 3, 5]

Output:
2: [[8.8] [-3.23333333333] [0.43333333333]]
1: 8.8+-3.23333333333*X+0.43333333333*X^2

Which means the polynomial fit for the given points is (approximately):

8.8 - 3.23333x + 0.43333x^2

(Screen shots of this example are shown below)

To get the length of the segment with the HP 50g, first note the xmin and xmax values, which are 1 and 6, respectively. This keystroke (in RPN Mode) should do the trick:

[ ' ] [ X ] [ENTER]
[RS] [COS] ( δ )
[LS] [ √ ] (x²)
1 [ + ] [ √ ]
1 [ENTER]
6 [ENTER]
3 [LS] [EVAL] (PRG) [F1] (STACK) [NXT] [F1] (ROLL)
[ ' ] [ X ] [ENTER]
[RS] [TAN] ( ∫ )

Length = 7.75662832829

My next blog entry will contain programs for the TI-84+ and Casio Prizm (fx series) calculators.

Have a great day everyone,

Eddie


This blog is property of Edward Shore. 2013

Length of a Polynomial Segment - Part 2 (TI-84+/Casio); From Coffee Klatch in San Dimas

Today I am coming to you from Coffee Klatch in San Dimas, CA. This is one of my most favorite places to hang out, and have great coffee and food.

Today's blog is the second part of the polynomial segment blog. Here is a link to Part 1.

---

To Recall:

A polynomial segment is a set of points that is connected by a polynomial.

To find a perfectly fit polynomial to fit n points:

1. Let X be the Vandermonde matrix of order (n-1) × (n-1) for the following:

2. Let y be the vector of y values.
3. The coefficient matrix is calculated as (X^T X)^-1 (X^T y). The result is a vector of coefficients of x^k where k = 0 to n-1.

On most graphing (and some scientific calculators), you can accomplish this task in two ways:

1. Use matrix operations. This is good for any size polynomial.
2. Take advantage of the curve fitting functions LinReg (to fit 2 points), QuadReg (to fit 3 points), CubicReg (to fit 4 points), and QuartReg (to fit 5 points). *Regression models and names vary for each calculator.

Once the polynomial, f(x) is found, find the arc length is found by:

The limits of the integral are the minimum and maximum x values, respectively.

---

Here are two programs that accomplish the task of finding a polynomial fit for 3, 4, or 5 points for the TI-84+ (and its family), and the Casio Prizm (and its fx-9xxx family).

---

TI-84+

POLYARC
1/12/2013
Fit a polynomial to a set of three, four, or five points.
Find the arc length of the polynomial from the end points.

Example 1: (0,0), (4,8), (8,6)
Results: y = -0.3125x^2 + 3.25
Length = 14.23920

Example 2: (0,3), (2,-1), (4,2), (6,0)
Results: y = -0.25x^3 + 2.375x^2 - 5.75x + 3
Length = 12.46834

Example 3: (0,0), (1,1), (2,0), (3,-1), (4,0.5)
Results: y ≈ -0.02083333x^4 + 0.20833333x^3 - 1.77083333x^2 + 2.54166666x
Length = 6.36142413

Program POLYARC (327 bytes)
Input "NO OF PTS(3-5)", N
If N<3 or N>5
1/0
N → dim(L1)
N → dim(L2)
For(K,1,N)
Disp K
Input "X=", X
Input "Y=", Y
X → L1(K)
Y → L2(K)
End
If N=3
Then
QuadReg L1, L2
fnInt(√(1+(2aX+b)²),X,min(L1),max(L1)) → L
Pause {a, b, c}
End
If N=4
Then
CubicReg L1, L2
Pause {a, b, c, d}
fnInt(√(1+(3aX² + 2bX + c)²), X, min(L1), max(L1)) → L
End
If N=5
Then
QuartReg L1, L2
Pause {a, b, c, d, e}
fnInt(√(1 + (4aX³ + 3bX² + 2cX + d)² , X, min(L1), max(L2)) → L
End
Disp "LENGTH="
Pause L

Where to find:
a, b, c, d, e: VARS, 5 (Statistics), right, right, [2, 3, 4, 5, and 6 respectively]
You need the lower case a, b, c, d, and e!

Note: No "bending backwards" is allowed

---

Casio Prizm

POLYARC
1/12/2013
Fit a polynomial to a set of three or four points.
Find the arc length of the polynomial from the end points.

Example 1: (0,0), (4,8), (8,6)
Results: y = -0.3125x^2 + 3.25
Length = 14.23920

Example 2: (0,3), (2,-1), (4,2), (6,0)
Results: y = -0.25x^3 + 2.375x^2 - 5.75x + 3
Length = 12.46834

Example 3: (0,0), (1,1), (2,0), (3,-1), (4,0.5)
Results: y ≈ -0.02083333x^4 + 0.20833333x^3 - 1.77083333x^2 + 2.54166666x
Length = 6.36142413

Program POLYARC (352 bytes):
"NO. OF POINTS(3,4,5)"? → N
N<3 Or N>5 ⇒ 1 ÷ 0
N → Dim List 1
N → Dim List 2
For 1 → K To N
K ◢
"X"? → List 1[K]
"Y"? → List 2[K]
Next
If N = 3
Then QuadReg List 1, List 2
∫( √(1 + (2aX + b)² ), Min(List 1), Max(List 1)) → L
{a, b, c} ◢
EndIf
If N=4
Then CubicReg List 1, List 2
∫( √(1 + (3aX ² + 2bX + c)² ), Min(List 1), Max(List 1)) → L
{a, b, c, d} ◢
EndIf
If N=5
Then QuartReg List 1, List 2
∫ ( √(1 + (4aX^3 + 3bX ² + 2cX + d)² ), Min(List 1), Max(List 2)) → L
{a, b, c, d, e} ◢
EndIf
"LENGTH="
L

How to Access:
a, b, c, d, and e: VARS, STAT (F3), GRAPH (F3), F1 (F2, F3, F4, and F5 respectively)
QuadReg: back to default menu, MENU, STAT, CALC, X^2
CubicReg: back to default menu, MENU, STAT, CALC, X^3
QuartReg: back to default menu, MENU, STAT, CALC, X^4

Note: No "bending backwards" allowed

---

Enjoy the day and I'll talk to you soon!

Eddie


This blog is property of Edward Shore. 2013

Arc length of sin(x) - Curve Approximation

Blog Entry #102

The Arc Length of a Sine Curve

The sine curve is one the most interesting curves in mathematics.

Let

y = a * sin x ,

Where a is the amplitude of the sine curve.

We can find the arc length of a curve between limits x1 and x2 by the integral:

x2
∫ √ (1 + (dy/dx)^2) dx
x1

For the sine curve:

y = a sin x

dy/dx = a cos x

(dy/dx)^2 = a^2 cos^x

And the arc length is:

x2
∫ √(1 + a^2 cos^x) dx
x1

There is no anti-derivative for √(1 + a^2 cos^x). Therefore, numerical methods must be used.

Finding an Approximate Curve

Using a TI nSpire CX CAS, I used the Spreadsheet, curve fitting, and graphing features to determine an approximate polynomial. The arc length is from the origin (0,0) to (π, 0).

Note: x1 = 0, x2 = π

Here is a shot summary of what I did:

1. Created a spreadsheet with the following columns:

Column A: A sequence of numbers from 0.25 to 5 in increments of 0.25. The resulting list is named amplist.

Column B: Use the nSpire's arcLen function to get the arc length of the sine curve using amplist as the values for a. This list is named arc1.

2. Pressing the menu key allowed me to access the Statistics menu. Using the Stat Calculations option, I used different types of regression, including power and quartic regression. What I was looking for was which regression had the best coefficient of determination (R^2). In general, the closer R^2 is to 1, the better the fit.

I ended up choosing the quartic regression (4th degree polynomial) with R^2 ≈ 0.9999919.

The approximate polynomial is

y = .0081196317102889 x^4 - .11577326164517 x^3 + .63914882375794 x^2 + .2071162669684 x + 3.0881429428239

This polynomial was save to the function f1(x).

3. I created a graphs page and made two plots:

* Scatter plot where x = amplist, y = arc1. (The dots in red)
* The function f1(x) (see step 2). (The curve in blue)

Setting zoom to fit the data, the curve looks like a good fit.

You can create a similar graph with the Data & Statstics module, but I thought I would be different this time.

How good of a fit is the polynomial?

4. I went back to the spreadsheet and added two more columns.

Column C: est1 = f1(amplist). (estimate arc lengths)
Column D: err1 = abs(arc1 - est1)

By scrolling down Column D, the quartic polynomial was accurate in estimating the arc length of the sine curve from 0 to π to at least two decimal places.

Conclusion

We have been looking to find the arc length of the curve y = a sin x from x = 0 to x = π.

The exact value is:

π
∫ √ (1 + a^2 cos^2 x ) dx
0

However, a good estimate can be found (to 2-3 decimal places) with the polynomial:

y = .0081196317102889 x^4 - .11577326164517 x^3 + .63914882375794 x^2 + .2071162669684 x + 3.0881429428239


Thanks as always and talk to you soon!


Eddie



This blog is property of Edward Shore. © 2012