Showing posts with label least squares. Show all posts
Showing posts with label least squares. Show all posts

Saturday, September 24, 2022

LSQ2: An update to LSQ (Casio fx-9750GIII, TI-84 Plus CE)

LSQ2:  An update to LSQ (Casio fx-9750GIII, TI-84 Plus CE)


Least Square Matrix and Correlation


The program LSQ2 fits the allows to fit data to a function with minimal error possible. 


Multiple Linear Regression:

f(x1, x2, x3, ...) = b0 + b1 * x1 + b2 * x2 + b3 * x3 + ....


Polynomial Regression:

f(x) = b0 + b1 * x + b2 * x^2 + b3 * x^3 + ...


General:

f(x) = b0 + b1 * g1(x) + b2 * g2(x) + b3 * g3(x) + ...


f(x1, x2, x3, ...) = b0 * g0(x1, x2, x3, ...) + b1 * g1(x1, x2, x3, ...) + ...


The Matrices X, Y, and B



X is your data matrix and is set up as columns:


[  g0(x),  g1(x),  g2(x), ... ]


Where the function g(x) represents functions applied to every data point x_i.  


Example 1:   f(x) = b0 + b1 * x 


The columns of the data matrix are set up as:

[ 1,  x ]


A column of ones set up solving for a constant.


Example 2:  f(x) = b0 + b1 * x + b2 * x^2


The columns of the data matrix are set up as:

[ 1, x, x^2 ]


Example 3:  f(x0, x1) = b0 + b1 * x1 + b2 * x2


The columns of the data matrix are set up as:

[ 1, x1, x2 ]  (note, not x squared in this case)



Y is the answer matrix, of size n rows and 1 column.  There are n data points. We are fitting the function to y_i.


B is the coefficient matrix, consisting of values b0, b1, b2, ....



Simply put, to find B using the least squares method given the data points:


B = (X^T X)^-1 X^T Y


X^T is the transpose matrix of X



How well does the function fit?  


We can predict y values by multiplying X by B.  


P' = X B 



Determining Coefficient of Correlation:


r^2 = SSreg ÷ SStot = [ B^T X^T Y - (O Y)^2 ÷ n ] ÷ [ Y^T Y - (O Y)^2 ÷ n ]

where O is a ones matrix [[ 1, 1, 1, 1, ... ]] of size 1 x n.  



Casio fx-9750GIII Program:  LSQ2


From the text file:  


'ProgramMode:RUN

ClrText

"2022_-_07_-_19 EWS"

"LEAST SQUARES"

"_Mat _X"?->Mat X

"_Mat _Y"?->Mat Y

Dim Mat Y->List 26

List 26[1]->N

(Trn Mat X*Mat X)^-1*Trn Mat X*Mat Y->Mat B

"_Mat _B:"Disps

Mat BDisps

{1,N}->Dim Mat O

Fill(1,Mat O)

Mat O*Mat Y->Mat S

Mat S*Mat S/N->Mat S

(Trn Mat B*Trn Mat X*Mat Y)-Mat S->Mat R

Mat R*(Trn Mat Y*Mat Y-Mat S)^-1->Mat R

"R_^<2>_:"Disps

Mat R



Listing:


ClrText

"2022-07-19 EWS"

"LEAST SQUARES"

"Mat X"? → Mat X

"Mat Y"? → Mat Y

Dim Mat Y → List 26

List 26[1] → N

(Trn Mat X × Mat X)^-1 × Trn Mat X × Mat Y → Mat B

"Mat B:" ⊿

Mat B ⊿

{1, N} → Dim Mat O

Fill(1, Mat O)

Mat O × Mat Y → Mat S

Mat S × Mat S ÷ N → Mat S

(Trn Mat B × Trn Mat X × Mat Y) - Mat S → Mat R

Mat R × (Trn Mat Y × Mat Y - Mat S)^-1 → Mat R

"R^2:" ⊿

Mat R


Matrices:

Mat X:  data matrix, X

Mat Y:  answer matrix, Y

Mat B: coefficient matrix, B

Mat O: ones matrix

Mat S:  used for calculation

Mat R:  correlation




TI-84 Plus CE Program:  LSQ2  (TI-Basic)



Listing:

"2022-07-19 EWS"
ClrHome
Disp "LEAST SQUARES"
Input "[X]? ", [J]
Input "[Y]? ", [I]
dim([I]) → L6
L6(1) → N
([J]^T [J])^-1 [J]^T [I] → [B]
Disp "[B]: "
Pause [B]
{1,N} → dim([H])
Fill(1,[H])
[H] [I] → [G]
[G] [G] * N^-1 → [G]
[B]^T [J]^T [I] - [G] → [A]
[A] * ([I]^T [I] - [G])^-1 → [A]
Disp "R^2: "
Disp [A]

List:
L6:  [ 2nd ] [ 6 ]

Matrices:
[J]:  data matrix, X
[I]:  answer matrix, Y
[B]: coefficient matrix, B
[H]: ones matrix
[G]:  used for calculation
[A]:  correlation


Examples

Example 1:

Equation: y = b0 + b1 * x1 + b2 * x2

X = [ [ 1, 1, 3 ] [ 1, 2, 4 ] [ 1, 5, 6 ] [ 1, 7, 3 ] [ 1, 7, 2 ] ]
Y = [ [ 0.86 ] [ 0.89 ] [ 0.95 ] [ 0.98 ] [ 0.96 ] ]

Coefficients:  [ [ b0 ] [ b1 ] [ b2 ] ]
B = [ [ 0.8257514451 ] [ 0.01836705202 ] [ 5.953757225E-3 ] ]

Correlation: [ [ 0.9875030926 ] ]

Example 2:

Equation: y = b0 + b1 * x + b2 * x^2

X = [ [ 1, 1, 1^2 ] [ 1, 2, 2^2 ] [ 1, 3, 3^2 ] [ 1, 4, 4^2 ] [ 1, 5, 5^2 ] [ 1, 6, 6^2 ] ]
Y = [ [ 1000 ] [ 1294 ] [ 1511 ] [ 1233 ] [ 1006 ] [ 879 ] ]

Coefficients:
B = [ [ 681.7 ] [ 435.2107143 ] [ -69.30357143 ] ]

Correlation: [ [ 0.8119609681 ] ]


Summary

Function to fit:   
f(x1, x2, x3 ... ) = b0 + b1 * g1(x1, x2, x3, ...) + b2 * g2(x1, x2, x3, ...) + ...

X = data matrix
Y = answer matrix, size n x 1
B = coefficient matrix

Determining the Coefficients:   B = (X^T X)^-1 X^T Y

Predicting Values:  P = X B

Determining Coefficient of Correlation:

r^2 = SSreg ÷ SStot = [ B^T X^T Y - (O Y)^2 ÷ n ] ÷ [ Y^T Y - (O Y)^2 ÷ n ]
where O is a ones matrix [[ 1, 1, 1, 1, ... ]] of size 1 x n.  

Source

Abdi, Hervè.  "Multiple Correlation Coefficient"  Program in Cognition and Neurosciences   https://personal.utdallas.edu/~herve/Abdi-MCC2007-pretty.pdf   Retrieved July 17, 2022.  


All original content copyright, © 2011-2022.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author. 

Sunday, November 29, 2020

HP Prime: Linear Exponential Combination Fit

 HP Prime:  Linear Exponential Combination Fit


The program LINEXPREG attempts to fit bivariate data to the curve:


y = a + b * x + c * e^x


The LSQ (least square) function is used.  The output is a list of three matrices:



*  A matrix of coefficients:  [ [ a ] [ b ] [ c ] ] 


*  A matrix of y values entered


*  A matrix of predicted y values



HP Prime Program LINEXPREG


EXPORT LINEXPREG(lx,ly)

BEGIN

// 2020-11-17 EWS

// x list, y list

LOCAL n,lx2,lx3,mx,my,k,mr,mq;

n:=SIZE(lx);

lx2:=e^(lx);

lx3:={};

FOR k FROM 1 TO n DO

lx3:=CONCAT(lx3,{1,lx(k),lx2(k)});

END;

mx:=list2mat(lx3,3);

my:=list2mat(ly,1);

mq:=LSQ(mx,my);

mr:=mq(1,1)+mq(2,1)*lx+

mq(3,1)*e^(lx);

mr:=list2mat(mr,1);

RETURN {mq,my,mr};

END;


Example





x list:  {0, 1, 2, 3, 4, 5}

y list:  {2, 7, 13, 18, 26, 34}


LINEXPREG({0, 1, 2, 3, 4, 5},{2, 7, 13, 18, 26, 34})


Results:  {coefficients, y values, predicted y values}


coefficients:

[ [ 1.74935499143 ]

[ 5.39455446221 ]

[ 3.66584105034E-2 ] ] 


y values:

[ [ 2 ]

[ 7 ]

[ 13 ]

[ 18 ]

[ 26 ]

[ 34 ] ]


predicted y values:

[ [  1.78601340193 ]

[ 7.24355734477 ]

[ 12.8093349675 ]

[ 18.6693222357 ]

[ 25.3290542368 ]

[ 34.1627178129 ] ] 


Equation:

y = 1.74935499143 + 5.39455446221 * x + 3.66584105034E-2 * e^x


On to the last month of 2020...


Eddie


All original content copyright, © 2011-2020.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author. 


Wednesday, April 12, 2017

Approximating y = sin x (0 to 90 degrees, 0 to pi/2 radians)



Approximating y = sin x (0 to 90 degrees, 0 to pi/2 radians)


Calculation

I used the HP Prime to fit polynomials to data by using the vandermonde and LSQ commands.

vandermonde(vector):  creates the matrix consisting of rows [ 1, (n_i), (n_i)^2, (n_i)^3, …, (n_i)^(n-1) ] where the vector has n elements for each i.  Available from the Math-Matrix-Create submenu.

LSQ(X, y):  returns the coefficients [ [a_0], [a_1], [a_2], … , [a_n] ], the minimum norm least squares vector from the system X*a=y.  In this case, the vector a can also be estimated by the operation (X^T X)^-1 X^T y.  Available from the Math-Matrix-Factorize submenu.

Case 1:  Four Points of Data

Fit cubic polynomial to the following data.  Remember we are working with the range of 0° ≤ x ≤ 90° or 0 ≤ x ≤ Ď€/2 radians

X
sin(x)
0
30°
1/2 = 0.5
60°
√3/2 ≈ 0.866025403785
90°
1

Example 1:  Full precision approximation

y = 1.78098409219E-2 x – 1.99435471445E-5 x^2 – 6.05408712068E-7 x^3

Note:  E-n stands for 10^(-n).

Absolute maximum error, data measured in 5° tick marks:  0.002371558429
Absolute maximum error, data measured in 1° tick marks:  0.002392465593



Example 2:  If we round all coefficients to six digits, we get:

y = 0.017810 x – 0.000020 x^2 – 0.000001 x^3

Absolute maximum error, data measured in 5° tick marks:  0.2881 

In comparison between Example 1 and Example 2, keeping the decimal places makes a big difference. 



Example 3:  Let’s see if we do better with using radians instead of degrees.

X
sin(x)
0
0
Ď€/6
1/2 = 0.5
Ď€/3
√3/2 ≈ 0.866025403785
Ď€/2
1

y = 1.02042871863 x – 0.065470803224 x^2 – 0.113871899065 x^3

Maximum absolute error, data measured in π/36 radian tick marks: 0.002371558428
Maximum absolute error, data measured in π/180 radian tick marks: 0.002392465661

The maximum error in Example 3 matched Example 1.

Example 4:  Round the coefficients to six digits:

y = 1.020429 x – 0.065471 x^2 – 0.113872 x^3

Maximum absolute error, data measured in Ď€/36 radian tick marks:  0.002371292924
(2 decimal places)

The difference between Example 2 and Example 4 is that less information is “lost” when rounding the coefficients to a set amount of digits, in this case, 6.

Case 2:  Seven Points of Data

x (degrees)
x (radians)
sin(x)
0
0
15°
Ď€/12
(√6 - √2)/4 ≈ 0.258819045102
30°
Ď€/6
1/2 = 0.5
45°
Ď€/4
√2/2 ≈ 0.707106781185
60°
Ď€/3
√3/2 ≈ 0.866025403785
75°
5Ď€/12
(√6 + √2)/4 ≈ 0.965925826288
90°
Ď€/2
1

Example 5:  Full Precision – Degrees

y = 1.74539026131E-2 x – 1.00636910048E-7 x^2 – 8.79821584089E-7 x^3 – 1.93973162632E-10 x^4 + 1.66950110845E-11 x^5 – 2.72879457224E-14 x^6

Absolute maximum error, data measured in 5° tick marks:  1.2072093E-6
(pretty darn good)

Example 6:  Full Precision – Radians

y = -6.98792887569E-12 + 1.0000349642 x – 3.30435720904E-4 x^2 – 0.165486304503 x^3 – 2.09062241321E-3 x^4 + 1.03087220879E-2 x^5  - 9.65423458482E-4 x^6

Maximum absolute error, data measured in Ď€/36 radian tick marks:  0.00049429473
(3 decimal places)

I think the degrees approximation wins.  But what if we round the radians version (Example 6) to 6 decimal places?

Example 7:  6 Decimal Places – Radians

y = 1.000035 x – 0.00033 x^2 – 0.165486 x^3 – 0.002091 x^4 + 0.010309 x^5 – 0.000965 x^6

Maximum absolute error, data measured in Ď€/36 radian tick marks:  9.03153E-6

Interesting result here, and a nice one at that!



Eddie

This blog is property of Edward Shore, 2017

Earth's Radius by Latitude

Earth's Radius by Latitude Introduction: Calculating the Earth’s Radius In quick, general calculations, we assume that the...