Polynomial Regression Model

Fabian Farelo, Autar Kaw, Jamie Trahan
University of South Florida
United States of America
kaw@eng.usf.edu

Introduction

Given n data points, ( (Subscript[x, 1]), (Subscript[y, 1])), ((Subscript[x, 2]), (Subscript[y, 2])), ..., ( (Subscript[x, n]), (Subscript[y, n])), use least squares method to regress the data to a (m^th)order polynomial.
y = (Subscript[a, 0])  +  (Subscript[a, 1]) x  +  (Subscript[a, 2]) (x^2) + ... + (Subscript[a, m]) (x^m),   m < n                                                                                                                     (1)

The residual at each data point is given by
(Subscript[E, i]) = (Subscript[y, i])-(Subscript[a, 0])-(Subscript[a, i]) (Subscript[x, i ]).... -(Subscript[a, m])((Subscript[x, i])^m)                                                                                                                                         (2)

The sum of the square of the residuals is given by
Sr = (∑_(i = 1)^n) ((Subscript[E, i])^2)												                    (3)

To find the constants of the polynomial regression model, we put the derivatives with respect to (Subscript[a, i]) to zero, that is,
\(TraditionalForm\`∂Sr/∂\((a\_0)\)\  = \ \(∑\+\(i = 1\)\%n\([2 \((\(y\_i - a\_0 - \(a\_i\) x\_i\  ... \) . \ \(-a\_m\) 
                      x\_\(\(i\)\(\ \)\)\^m)\)]\) \((\(-1\))\)\  = \ 0\)\)							       (4.a)
\(TraditionalForm\`∂Sr/∂\((a\_1)\)\  = \ \(∑\+\(i = 1\)\%n\([2 \((\(y\_i - a\_0 - \(a\_i\) x\_i\  ... \) . \ \(-a\_m\) 
                      x\_\(\(i\)\(\ \)\)\^m)\)]\) \((\(-x\_i\))\)\  = \ 
          0\)\)							       (4.b)
.		.		.		.		.		.		.		.		.
.		.		.		.		.		.		.		.		.
\(TraditionalForm\`∂Sr/∂\((a\_m)\)\  = \ \(∑\+\(i = 1\)\%n\([2 \((\(y\_i - a\_0 - \(a\_i\) x\_i\  ... \) . \ \(-a\_m\) 
                      x\_\(\(i\)\(\ \)\)\^m)\)]\) \((\(-x\_i\^m\))\)\  = \ 
          0\)\)						                  (4.m)

Setting those equations in matrix form gives
(({{n, ∑(Subscript[x, i]), ..., ∑((Subscript[x, i])^m)}, {∑(Subscript[x, i]), ∑((Subscript[x, i])^2), ..., ∑((Subscript[x, i])^(m + 1))}, {..., ..., ..., ...}, {∑((Subscript[x, i])^m), ∑((Subscript[x, i])^(m + 1)), ..., ∑((Subscript[x, i])^(2  m))}}))(({{(Subscript[a, 0])}, {(Subscript[a, 1])}, {..}, {(Subscript[a, m])}}))=(({{∑(Subscript[y, i])}, {∑(Subscript[x, i]) (Subscript[y, i])}, {...}, {∑((Subscript[x, i])^m) (Subscript[y, i])}}))

The above simultaneous linear equations are solved for the (m +1) constants, (Subscript[a, 0]),  
      (Subscript[a, 1]),  ... (Subscript[a, m]). To learn more about polynomial regression see the Nonlinear Regression models worksheet.

Section 1: Input Data

Below are the input parameters to begin the simulation. This is the only section that requires user input. Once the values are entered, Mathematica will generate a polynomial regression model for the given data set.

Number of Equations, n

n = 11

Array of x values, "X"

X = {80, 40, 
        0, -40, -80, -120, -160, -200, -240, -280, -320}

Array of y values, "Y"

Y = {6.47, 6.24, 6. , 5.72, 5.43, 5.09, 4.72, 4.30, 3.83, 3.33, 
        2.76}

Desired order of polynomial regression model

OrderPoly = 1

Enter the lowest order of polynomial to check for optimum order:

LowOrder = 1

Enter the highest order of polynomial to check for optimum order.   NOTE: HighOrder must be less than or equal to n-1.

HighOrder = 6

Section 2: Polynomial Regression Model

In this section, the coefficient matrix "M" and right hand side vector "B" are calculated and subsequently used to determine the solution vector that contains the coefficients of the polynomial model  (Subscript[a, 0]),  
      (Subscript[a, 1]),  ... (Subscript[a, m]).

Determining the coefficient matrix, "M"

Creating the matrix size according to the order of the polynomial:

M = Array[0, {OrderPoly + 1, OrderPoly + 1}];

Determining each value of the first row of matrix "M":

M[[1, 1]] = n;
For[i = 2, i ≤ OrderPoly + 1, i++, 
        M[[1, i]] = 0; 
For[j = 1, 
          j ≤ n, j++, 
M[[1, i]] = 
            N[M[[1, i]] + X[[j]]^(i - 1)]]];

Calculating the remaining values of the coefficient matrix:

For[i = 1, 
      i ≤ OrderPoly + 1, i++, 
For[k = 2, 
        k ≤ OrderPoly + 1, k++, 
        M[[k, i]] = 0; 
For[j = 1, 
          j ≤ n, j++, 
M[[k, i]] = 
            N[M[[k, i]] + X[[j]]^(i + k - 2)]]]]

We now have the coefficient matrix, "M":

M // MatrixForm

Determining the right hand side vector, "B"

Creating the vector size according to the order of the polynomial model.

B = Array[0, {OrderPoly + 1}];

Finding the value of the array

B[[1]] = 0;

For[i = 1, 
        i ≤ n, i++, 
B[[1]] = 
          B[[1]] + Y[[i]]];
For[i = 2, i ≤ OrderPoly + 1, i++, 
      B[[i]] = 0; 
For[j = 1, 
        j ≤ n, j++, 
B[[i]] = 
          N[B[[i]] + (X[[j]]^(i - 1))*
                Y[[j]]]]]

Now, the right hand side vector "B"

B // MatrixForm

Calculating the coefficients of the model

Merging the matrices "M" and "B" to use Mathematica to reduce it:

F = Transpose[Append[Transpose[M], B]];

Selecting the last column with the coefficient's values.
This array contains the coefficients of the model as a = [(Subscript[a, 0]),  (Subscript[a, 1]),   ...  ,  (Subscript[a, m])].

a = RowReduce[F];
a = a[[All, OrderPoly + 2]]

The polynomial regression model is as follows:

For[i = 1, 
        i ≤ n, i++, 
y = 
          a[[1]]; 
For[j = 1, 
          j ≤ OrderPoly, j++, 
y = 
            y + a[[(j + 1)]]*x^j]];
Print["y= " y]

The plot shows the data points as well as the regression model.

Observed = 
        Table[{X[[i]], Y[[i]]}, {i, 1, 
            n}];
predicted = y;
ttl = 
      Print["Polynomial Regression Model of order ", 
        OrderPoly]
pt = 
        ListPlot[Observed, PlotStyle -> PointSize[0.01], 
          DisplayFunction -> Identity];
lin = 
        Plot[y, {x, Min[X], Max[X]}, 
          DisplayFunction -> Identity];
\(Show[pt, lin, 
      AxesLabel -> {"X", \*""y=(Subscript[a, 0])+(Subscript[a, 1])x+(Subscript[a, 2])(x^2)+...(Subscript[a, m])(x^m)""}, DisplayFunction -> $DisplayFunction]\)

Section 3: Optimum Order

In this section we will determine the optimum order of the polynomial model by plotting the variance defined as 

((Subscript[S, r])/(n - (m + 1)))

as a function of m, where n is the number of data points, (Subscript[S, r]) is the sum of the square of residuals and m is the order of the polynomial. The optimum order is considered as to be the one where the value of the variance ((Subscript[S, r])/(n - (m + 1))) is minimum or where its value is significantly decreasing. 
In the following procedure, a polynomial regression model is calculated for each order specified in the LowOrder to HighOrder range. Mathematica then computes the variance of each model.The worksheet does not choose the order of the optimum polynomial for regression for you.  Look at the plot of the variance as a function of the order of the polynomial.  The optimum polynomial is one after which there is no statistical significant decrease in the variance.  

 Many a times, the variance may show signs of decreasing and then increasing as a function of the order of the polynomial regression model.  Such increases in the variance are normal as the variance is calculated as the ratio between the sum of the squares of the residuals and the difference between the number of data points and number of constants of the polynomial model.  Both the numerator and denominator decrease as the order of the polynomial is increased.  However, as the order of the polynomial increases, the coefficient matrix in the calculation of the constants of the model becomes more ill-conditioned. This ill-conditioning of the coefficient matrix results in fewer significant digits that can be trusted to be correct in the coefficients of the polynomial model, and hence artificially amplify the value of the variance.

Sr = Array[0, {HighOrder}];
  

For[m = LowOrder, 
      m ≤ HighOrder, m++, 
  
M = 
        Array[0, {m + 1, m + 1}]; 
M[[1, 1]] = 
        n; 
For[i = 2, i ≤ m + 1, i++, 
        M[[1, i]] = 0; 
For[j = 1, 
          j ≤ n, j++, 
          M[[1, i]] = 
            N[M[[1, i]] + 
                X[[j]]^(i - 1)]]]; 
For[i = 1, 
        i ≤ m + 1, i++, 
For[k = 2, 
          k ≤ m + 1, k++, 
          M[[k, i]] = 0; 
For[j = 1, 
            j ≤ n, j++, 
M[[k, i]] = 
              N[M[[k, i]] + 
                  X[[j]]^(i + k - 2)]]]]; 
  
q = 
        Array[0, {m + 1}]; 
q[[1]] = 
        0; 
For[i = 1, 
        i ≤ n, i++, 
q[[1]] = 
          q[[1]] + Y[[i]]]; 
For[i = 2, 
        i ≤ m + 1, i++, 
        q[[i]] = 0; 
For[j = 1, 
          j ≤ n, j++, 
q[[i]] = 
            N[q[[i]] + (X[[j]]^(i - 1))*
                  Y[[j]]]]]; 
  
F = 
        Transpose[Append[Transpose[M], q]]; 
a = 
        RowReduce[F]; 
a = 
        a[[All, m + 2]]; 
  
Sr[[m]] = 
        0; 
For[i = 1, i ≤ n, i++, 
        summ = 0; 
For[j = 1, j ≤ m, j++, 
          summ = summ + 
              a[[j + 1]]*X[[i]]^j]; 
Sr[[m]] = 
          Sr[[m]] + (Y[[i]] - (a[[1]] + 
                      summ))^2];]

      
var = Array[0, {HighOrder}];
For[i = LowOrder, i ≤ HighOrder, i++, 
        var[[i]] = 0; 
var[[i]] = 
          N[Sr[[i]]/(n - (i + 1))]];
Print["var =",  var];

\(ListPlot[var, PlotStyle -> PointSize[0.02], 
      PlotJoined -> True, 
      PlotLabel -> "Optimum Order of Polynomial", 
      AxesLabel -> {"Order of Polynomial, m", "Sr/[n-(m+1)]"}]\)

Conclusion

Using Mathematica, we are able to regress a given data set to a polynomial model of the (m^th) order.

Question 1: Water is flowing through a pipe of radius 0.5 ft and flow velocity, v measurements are made from the center of the wall of the pipe as follows:

({{Radial Location,
r (ft), Velocity, 
v (ft/s)}, {0, 10}, {0.08, 9.7}, {0.17, 8.9}, {0.25, 7.5}, {0.32, 5.6}, {0.42, 3.1}, {0.50, 0}})

a) Regress the data to 
    v = (Subscript[a, 0])(1 - ((r^2)/(a^2)))
    where a is the radius of the pipe.
b) Find the flow rate, \(TraditionalForm\`Q\& . \) through the pipe. (Hint: \(TraditionalForm\`Q\& . \) = ∫2πrvdr, r=0..a)

Question 2: Thermal expansion coefficient of steel varies with temperature as given in the table below.

({{          Temperature, T          , Thermal expansion coefficient, α*E - 06}, {80, 6.47}, {60, 6.36}, {40, 6.24}, {20, 6.12}, {0, 6.00}, {-20, 5.86}, {-40, 5.72}, {-60, 5.58}, {-80, 5.43}, {-100, 5.28}, {-120, 5.09}, {-140, 4.91}, {-160, 4.72}, {-180, 4.52}, {-200, 4.30}, {-220, 4.08}, {-240, 3.83}, {-260, 3.58}, {-280, 3.33}, {-300, 3.07}, {-320, 2.76}, {-340, 2.45}})

a) Regress the data to a second order polynomial.
   α = (Subscript[a, 0])  +  (Subscript[a, 1]) T  +  (Subscript[a, 2]) (T^2)
b) Find the optimum order of polynomial for the regression model. 
c) Find the reduction in the diameter of a steel cylinder of diameter 12.5" if it is cooled from a room temperature of 80°F to -180°F.

References

[1]Autar Kaw, Holistic Numerical Methods Institute, http://numericalmethods.eng.usf.edu/nbm, See
How does Nonlinear Regression work?