Using Matrix Equations to Find Regression Polynomials

Date: 06/28/2006 at 10:00:16
From: mary ann
Subject: How do I find a Cubic Fit

I would like to know if there is any formula for a cubic regression. 
Generally in the books that I have looked at, they say make a scatter
plot of the data.  Is there some definite formula that I can use for
determining the coefficients of a cubic regression model?

Date: 06/28/2006 at 12:23:40
From: Doctor Fenton
Subject: Re: How do I find a Cubic Fit

Hi Mary Ann,

Thanks for writing to Dr. Math.  Yes, there is a formula.  If you have
never seen any value to the concept of matrices, this subject should
quickly convince you.  With matrix notation and concepts, it is easy
to write a (matrix) equation for the coefficients of the regression
polynomial, which can be solved numerically (it usually isn't worth
solving analytically).

To fit a cubic y = a_0 + a_1*x + a_2*x^2 + a_3*x^3

to a list of points (x_1,y_1),(x_2,y_2),...,(x_n,y_n)  with
n >=4,  you can write the equations

    a_0 + a_1*x_1 + a_2*(x_1)^2 + a_3*(x_1)^3 = y_1
    a_0 + a_1*x_2 + a_2*(x_2)^2 + a_3*(x_2)^3 = y_2
    a_0 + a_1*x_3 + a_2*(x_3)^2 + a_3*(x_3)^3 = y_3
        :         :             :                :
    a_0 + a_1*x_n + a_2*(x_n)^2 + a_3*(x_n)^3 = y_n

which can be written as a matrix equation

    AX = b

with     [ 1 x_1 (x_1)^2 (x_1)^3 ]       [ a_0 ]        [ y_1 ]
         [ 1 x_2 (x_2)^2 (x_2)^3 ]       [ a_1 ]        [ y_2 ]
     A = [ 1 x_3 (x_3)^2 (x_3)^3 ]   x = [ a_2 ] ,  b = [ y_3 ]
            :    :      :      :         [ a_3 ]           :
         [ 1 x_n (x_n)^2 (x_n)^3 ] ,                    [ y_n ]  .

You can show that if this system of n equations in 4 unknowns has rank 
4 (at least 4 linearly independent rows), then the system of "normal 
equations" (as they are called)

    A^t A X = A^t b 

where A^t is the transpose of A, so that A^tA is a 4x4 matrix, gives
the solution of the regression equation.  

For example, in the linear case (fitting a regression line 
y = a_0 + a_1*x), the normal equations are

 [  1   1   1  ...  1  ][ 1 x_1 ] [ a_0 ]   [  1   1  ...  1  ][y_1]  
 [ x_1 x_2 x_3 ... x_n ][ 1 x_2 ] [ a_1 ] = [ x_1 x_2 ... x_n ][y_2]
                        [ 1 x_3 ]                              [y_3]
                           :  :                                  :
                        [ 1 x_n ]                              [y_n]

which becomes

               n                      n
   [ n        SUM(x_k)    ]        [ SUM y_k     ]
   [          k=1         ][a_0]   [ k=1         ]
   [  n        n          ][a_1] = [             ]
   [ SUM x_k  SUM (x_k)^2 ]        [  n          ]
   [ k=1      k=1         ]        [ SUM x_k*y_k ]    .
                                     k=1

From this, using Cramer's Rule, you can obtain the linear regression
formulas.

You can use Cramer's Rule, for example, to write out the explicit
formulas for the coefficients a_i, i=0,1,2,3, but they will be quite
complicated, since that would require the expansion of 4x4
determinants.  In practice, one either sets up the normal equations
and solves numerically, or uses an alternative formulation based on
matrix decompositions such as the QR or Singular Value Decomposition.

This approach works for any degree polynomial, so long as the rank of
A exceeds the degree of the polynomial plus 1.

If you have any questions, please write back and I will try to 
explain further.

- Doctor Fenton, The Math Forum
  http://mathforum.org/dr.math/