Nonlinear Regression with Data Linearization

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

ClearAll;

Introduction

This worksheet illustrates finding the constants of nonlinear regression models with data linearization. Three common nonlinear models are illustrated - 1) Exponential  2) Power  3) Saturation Growth.

Given n data points ( (Subscript[x, 1]), (Subscript[y, 1])), ( (Subscript[x, 2]), (Subscript[y, 2])), ( (Subscript[x, 3]), (Subscript[y, 3])), .... , ( (Subscript[x, n]),  (Subscript[y, n])), best fit one of the following models to the data.

Exponential: y = (ae^bx)
Power: y = (ax^b)
Saturation:  y = (ax/(b + x))

a and b are the constants of the regression model.

To learn more about Nonlinear Regression with data linearization, see the Nonlinear Regression worksheet.

Section 1: Input Data

Below are the input parameters to begin the simulation. This is the only section that requires user input. The user can change the values below and Mathematica will calculate the nonlinear regression model using data linearization for the corresponding model that the user chooses. 

Input Parameters:

Number of data points, n

n = 5

Array of x values, x

x = {10, 16, 25, 40, 60}

Array of y values, y

y = {94, 118, 147, 180, 230}

For exponential model call model type to be "exponential"
For power model assign the model type variable as "power"
For saturation growth model, assign the model type variable to be "growth"

modeltype = "exponential"

Section 2: Nonlinear Model with Data Linearization Procedures

The following manipulation must first be made to each model in order to linearize the data:

Exponential: ln(y) = ln(a) + bx       (taking the natural log of both sides)                                   (2.1)
Power: ln(y) = ln(a) + bln(x)           (taking the natural log of both sides)                                   (2.2)
Growth: (1/y)= (b/a) (1/x) + (1/a)                (taking the reciprocal of both sides then rearranging)          (2.3)

Once the data is linearized, substitutions are made so as to apply a direct solution approach using least squares regression method for a straight line. In this section, a procedure is written for each model. Once a specific nonlinear model is called for with the modeltype variable, the proper procedure will calculate the coefficients for that model using data linearization.
Note: See Nonlinear Regression notes for explanation of the model.

Exponential model procedure:

After data linearization, the following substitutions are made to Equation 2.1
 let   z = ln(y)
     (Subscript[a, 0]) = ln(a),  implying a = (Subscript[(e^a), 0])
      (Subscript[a, 1])= b
The data z versus x now takes the form of a linear model:
      
   z = (Subscript[a, 0])  + (Subscript[a, 1]) x                                                                                           (2.4)
   
In the exponential model procedure, least squares linear regression method is used to solve for the (Subscript[a, 0])and (Subscript[a, 1])coefficients which are then used to determine the original constants of the exponential model, a and b, where y = (ae^bx).

Off[General::spell]

ExponentialModel[x_, y_, n_] := 
        Module[{z, sumz, sumx, sumxz, sumxx, a1, a0}, 
z = 
            Array[0, {1, n}]; sumz = 0; sumx = 0; sumxz = 0; 
          sumxx = 0; 
For[i = 1, i ≤ n, i++, 
            z[[1, i]] = 
              N[Log[y[[i]]]]; 
      sumz = 
              sumz + z[[1, i]]; 
      sumx = 
              sumx + x[[i]]; 
      sumxz = 
              sumxz + x[[i]]*
                  z[[1, i]]; 
      sumxx = 
              sumxx + (x[[i]])^2;]; 
a1 = 
            N[((n*sumxz - sumx*sumz)/((n*sumxx) - sumx^2))]; 
a0 = N[(sumz/n) - a1*(sumx/n)]; 
a = 
            E^a0; 
b = a1; 
a; 
          b];

Power model procedure:

After data linearization, the following substitutions are made to Equation 2.2
let    z = ln(y)
      w = ln(x)
      (Subscript[a, 0]) =  ln(a),   implying a = (e^a0)
      (Subscript[a, 1])= b
 The data z versus w now takes the form of a linear model:
     
     z = (Subscript[a, 0])+ (Subscript[a, 1]) w                                                                                  (2.5)
 
In the power model procedure, least squares linear regression method is used to solve for the (Subscript[a, 0]) and (Subscript[a, 1]) coefficients which are then used to determine the originial constants of the power model, a and b, where y = (ax^b).

PowerModel[x_, y_, n_] := 
        Module[{z, w, sumz, sumw, sumwz, sumww}, 
z = 
            Array[0, {1, n}]; w = Array[0, {1, n}]; sumz = 0; sumw = 0; 
          sumwz = 0; 
          sumww = 0; 
For[i = 1, i ≤ n, i++, 
            z[[1, i]] = 
              N[Log[y[[i]]]]; 
w[[1, i]] = 
              N[Log[x[[i]]]]; 
sumw = 
              sumw + w[[1, i]]; 
sumz = 
              sumz + z[[1, i]]; 
sumwz = 
              sumwz + w[[1, i]]*z[[1, i]]; 
sumww = 
              sumww + (w[[1, i]])^2;]; 
a1 = 
            N[((n*sumwz - sumw*sumz)/(n*sumww - sumw^2))]; 
a0 = N[(sumz/n) - a1*(sumw/n)]; 
a = 
            E^a0; 
b = a1] ;

Saturation Growth model procedure:

After data linearization, the following substitutions are made to Equation 2.3
let    z = 1/y
       q = 1/x
     (Subscript[a, 0])  =  1/a,    implying a = 1/(Subscript[a, 0])
      (Subscript[a, 1])= b/a,      implying b = (Subscript[a, 1])/(Subscript[a, 0])
The data z versus q now takes the form of a linear model:

z = (Subscript[a, 0]) + (Subscript[a, 1]) q                                                                                     (2.6)

In the saturation growth model procedure, least squares linear regression method is used to solve for the (Subscript[a, 0]) and (Subscript[a, 1])coefficients which are then used to determine the original constants of the growth model, a and b, where y = (ax/(b + x)).

GrowthModel[x_, y_, n_] := 
        Module[{z, q, sumq, sumz, sumqz, sumqq, a1, 
            a0}, 
z = Array[0, {1, n}]; 
          q = Array[0, {1, n}]; sumq = 0; sumz = 0; sumqz = 0; 
          sumqq = 0; 
For[i = 1, i ≤ n, i++, 
            z[[1, i]] = 
              N[1/y[[i]]]; 
q[[1, i]] = 
              N[1/x[[i]]]; 
sumq = 
              sumq + q[[1, i]]; 
sumz = 
              sumz + z[[1, i]]; 
sumqz = 
              sumqz + q[[1, i]]*z[[1, i]]; 
sumqq = 
              sumqq + (q[[1, i]])^2;]; 
a1 = 
            N[((n*sumqz - sumq*sumz)/(n*sumqq - sumq^2))]; 
a0 = N[(sumz/n) - a1*(sumq/n)]; 
a = 
            N[1/a0]; 
b = N[a1/a0];];

Section 3: Finding the constants of the model

The proper regression model is assigned below, returning the constants of the model that was chosen with the modeltype variable.

If[modeltype == "exponential", ExponentialModel[x, y, n]; 
        f[x_] := N[a*E^(b*x)]];
If[modeltype == "power", PowerModel[x, y, n]; 
        f[x_] := a*x^b];
If[modeltype == "growth", GrowthModel[x, y, n]; 
        f[x_] := (a*x)/(b + x)];
Print["The constants of the ", modeltype, " model are a = ", 
      a, " and b = ", b]

Below, the nonlinear model is plotted versus the data set.

Datatable = 
        Table[{x[[i]], y[[i]]}, {i, 1, 
            n}];
title = 
        Print[modeltype, " Regression Model, y vs x"];
points = 
        ListPlot[Datatable, PlotStyle -> PointSize[0.02], 
          DisplayFunction -> Identity];
lin = 
        Plot[f[X], {X, Min[x], Max[x]}, 
          DisplayFunction -> Identity];
Show[points, lin, 
        DisplayFunction -> $DisplayFunction];

Conclusion

Mathematica helped us apply our knowledge of linear regression method to regress a given data set to a nonlinear model

Question 1: You are working for Valdez SpillProof Oil Company as a petroleum engineer.  Your boss is asking you to estimate the future life of an oil well.  The analysis used in the industry is called the decline curve analysis where the barrels of oil produced per unit time are plotted against time, and the curve is extrapolated.  One of the standard curves used is harmonic decline model, that is
				q =     b   
				          1 + at
where q is the rate of production and t is the time, b and a are the constants of the regression model. 
({{time (t),  Month, 1, 5, 11, 16, 18}, {\(Rate\ of\ Production\ \((q)\), 
Barrels\ per\ Day\), 258, 185, 109, 84, 71}})

a)  Find the constants of the regression model. Hint: You are allowed to linearize the data if possible.
b)  Find the total life of an oil field if 5 barrels per day is considered the production at  which the field needs to be abandoned for further production.
c) What does b stand for?

Question 2: It is desired to obtain a functional relationship between the mass density  r of air and the altitude h above the sea level for the dynamic analysis of bodies moving within earth's atmosphere.  Use the approximation r = (Subscript[k, 1]) (e^(-(Subscript[k, 2]) h)) to fit the data given below by regression analysis.  Find the constants k1 and k2.  You are allowed to linearize the data.
({{Altitude

            (kilometers), Mass Density

            (kg/(m^3))}, {0.32, 1.15}, {0.64, 1.10}, {1.28, 1.05}, {1.60, 0.95}})

Question 3: It is suspected from theoretical considerations that the rate of flow from a fire hose is proportional to some power of the nozzle pressure.  Determine whether the speculation is true.  What is the exponent of the data.  Assume pressure data is more accurate.  You are allowed to linearize the data.

({{Flow rate (gallons/min) F, 94, 118, 147, 180, 230}, {Pressure (psi) p, 10, 16, 25, 40, 60}})

References

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