Chapter 4

Curve Fitting

In this section, we demonstrate the ideas of fitting curves to data using the X-Files, Cost of Advertising, Black Bear, and Corn Storage data.

Mathematica Background

This Chapter relies on two commands, the Fit command and the Regress command.  The Fit command fits any intrinsically linear function to data.  Its syntax is Fit[data, {basis functions}, independent variables].  If an ANOVA or regression table is desired one must use the Regress command instead.  The Regress command is in the Statistics`LinearRegression package.

Reopening the X Files

The following data set  is the number of viewers (in millions) for the first time airing of the television show The X Files for each week of its first  season, second season, and both seasons combined. The data is week number versus the number of the viewers (in millions).

season1
       = {{1, 7.4}, {2, 6.9}, {3, 6.8}, {4, 5.9}, {5, 6.2}, {7, 5.6}, {8, 
          5.6}, {9, 6.2}, {10, 6.1}, {11, 5.1}, {14, 6.4}, {15, 6.4}, {18, 
          6.2}, {20, 6.8}, {22, 7.2}, {23, 6.8}, {23, 5.8}, {24, 7.1}, {28, 
          7.2}, {30, 7.5}, {32, 8.1}, {33, 7.7}, {35, 7.4}, {36, 8.3}};

season2
       = {{53, 9.8}, {54, 9.3}, {55, 8.7}, {56, 8.2}, {57, 8.5}, {58, 9.2}, {
          60, 9}, {61, 9.1}, {62, 8.6}, {65, 9.9}, {66, 8.5}, {69, 9.7}, {70, 
          8.8}, {72, 10.2}, {73, 10.8}, {74, 9.8}, {75, 10.7}, {76, 9.6}, {
          78, 10.2}, {80, 9.8}, {82, 7.9}, {84, 8.5}, {85, 8.1}, {86, 9}, {
          87, 9.6}};

seasons
       = {{1, 7.4}, {2, 6.9}, {3, 6.8}, {4, 5.9}, {5, 6.2}, {7, 5.6}, {8, 
          5.6}, {9, 6.2}, {10, 6.1}, {11, 5.1}, {14, 6.4}, {15, 6.4}, {18, 
          6.2}, {20, 6.8}, {22, 7.2}, {23, 6.8}, {23, 5.8}, {24, 7.1}, {28, 
          7.2}, {30, 7.5}, {32, 8.1}, {33, 7.7}, {35, 7.4}, {36, 8.3}, {53, 
          9.8}, {54, 9.3}, {55, 8.7}, {56, 8.2}, {57, 8.5}, {58, 9.2}, {60, 
          9}, {61, 9.1}, {62, 8.6}, {65, 9.9}, {66, 8.5}, {69, 9.7}, {70, 
          8.8}, {72, 10.2}, {73, 10.8}, {74, 9.8}, {75, 10.7}, {76, 9.6}, {
          78, 10.2}, {80, 9.8}, {82, 7.9}, {84, 8.5}, {85, 8.1}, {86, 9}, {
          87, 9.6}};

We begin by creating a scatterplot of the data from the first season.

season1Plot = ListPlot[season1,
AxesLabel->{"Week Number","Viewers in Millions"},
PlotStyle->{PointSize[.02]}];

Next we fit a line to the first season data and overlay the line and th edata.

f[x_]=Fit[season1,{1,x},x]

season1ext = Plot[f[x],{x,0,36},
AxesLabel->{"Week Number","Viewers in Millions"},DisplayFunction->Identity];
Show[season1Plot,season1ext,DisplayFunction->$DisplayFunction];

Finally we compute (R^2) directly from the definition.

data = season1;

SSres = 
      (∑_(n = 1)^(Length[data]))(((
            data[[n
                  ]][[2
                ]] - 
              f[data[[n
                    ]][[1
                  ]]]))^2)

mean = 
      (((∑_(n = 1)^(Length[data]))data[[n
              ]][[2
            ]])/Length)[data]

SSreg = 
      (∑_(n = 1)^(Length[data]))(((
            f[data[[n
                    ]][[1
                  ]]] - mean))^2)

SStot = SSres + SSreg

(SSreg/SStot)

Thus R\(TraditionalForm\`\^2\) is .41, so this linear model accounts for 41% of the variance in the data.

Next we will see how the first season model does at explaining the variation seen in both seasons.

seasonsPlot = ListPlot[seasons,
AxesLabel->{"Week Number","Viewers in Millions"},
PlotStyle->{PointSize[.02]}];

seasonsext = Plot[f[x],{x,0,90},
AxesLabel->{"Week Number","Viewers in Millions"},DisplayFunction->Identity];
Show[seasonsPlot,seasonsext,DisplayFunction->$DisplayFunction];

data = seasons;

SSres = 
      (∑_(n = 1)^(Length[data]))(((
            data[[n
                  ]][[2
                ]] - 
              f[data[[n
                    ]][[1
                  ]]]))^2)

mean = 
      (((∑_(n = 1)^(Length[data]))data[[n
              ]][[2
            ]])/Length)[data]

SSreg = 
      (∑_(n = 1)^(Length[data]))(((
            f[data[[n
                    ]][[1
                  ]]] - mean))^2)

SStot = SSres + SSreg

(SSreg/SStot)

The first season model explains 75% of the variation in the data for both seasons.  Interestingly, this is better than the amount of variation in the first season's data explained by the first season's model.

Now let's try finding R^2 of the first season's model with the data of the second season:

data = season2;

SSres = 
      (∑_(n = 1)^(Length[data]))(((
            data[[n
                  ]][[2
                ]] - 
              f[data[[n
                    ]][[1
                  ]]]))^2)

mean = 
      (((∑_(n = 1)^(Length[data]))data[[n
              ]][[2
            ]])/Length)[data]

SSreg = 
      (∑_(n = 1)^(Length[data]))(((
            f[data[[n
                    ]][[1
                  ]]] - mean))^2)

SStot = SSres + SSreg

(SSreg/SStot)

The first season model explains 44% of the variation in the data for both seasons.  Interestingly, this is better than the amount of variation in the first season's data explained by the first season's model.

Let's try this again:

We begin by creating a scatterplot of the data from the first season.

season2Plot = ListPlot[season2,
AxesLabel->{"Week Number","Viewers in Millions"},
PlotStyle->{PointSize[.02]}];

Next we fit a line to the second season data and overlay the line and th edata.

f2[x_]=Fit[season2,{1,x},x]

season2ext = Plot[f2[x],{x,53,87},
AxesLabel->{"Week Number","Viewers in Millions"},DisplayFunction->Identity];
Show[season2Plot,season2ext,DisplayFunction->$DisplayFunction];

Finally we compute (R^2) directly from the definition.

data = season2;

SSres = 
      (∑_(n = 1)^(Length[data]))(((
            data[[n
                  ]][[2
                ]] - 
              f2[data[[n
                    ]][[1
                  ]]]))^2)

mean = 
      (((∑_(n = 1)^(Length[data]))data[[n
              ]][[2
            ]])/Length)[data]

SSreg = 
      (∑_(n = 1)^(Length[data]))(((
            f2[data[[n
                    ]][[1
                  ]]] - mean))^2)

SStot = SSres + SSreg

(SSreg/SStot)

Thus R\(TraditionalForm\`\^2\) is 0.00981683, so this linear model accounts for 0.981683% of the variance in the data.  That's almost 1% of the variation! Yuck!

Next we will see how the second season model does at explaining the variation seen in both seasons.

seasonsPlot = ListPlot[seasons,
AxesLabel->{"Week Number","Viewers in Millions"},
PlotStyle->{PointSize[.02]},DisplayFunction->Identity];
seasonsext = Plot[f2[x],{x,0,90},DisplayFunction->Identity,
AxesLabel->{"Week Number","Viewers in Millions"}];
Show[seasonsPlot,seasonsext,DisplayFunction->$DisplayFunction];

data = seasons;

SSres = 
      (∑_(n = 1)^(Length[data]))(((
            data[[n
                  ]][[2
                ]] - 
              f2[data[[n
                    ]][[1
                  ]]]))^2)

mean = 
      (((∑_(n = 1)^(Length[data]))data[[n
              ]][[2
            ]])/Length)[data]

SSreg = 
      (∑_(n = 1)^(Length[data]))(((
            f2[data[[n
                    ]][[1
                  ]]] - mean))^2)

SStot = SSres + SSreg

(SSreg/SStot)

The first season model explains 29% of the variation in the data for both seasons.  Interestingly, this is better than the amount of variation in the second season's data explained by the second season's model.

Cost of Advertising

This section shows the complete analysis used in the Cost of Advertising example.

The annual advertising expenditures (in millions) in the U.S. from 1970 to 1989.

CostOfAdv
       = {{70, 19550}, {71, 20700}, {72, 23210}, {73, 24980}, {74, 26620}, {
          75, 27900}, {76, 33300}, {77, 37440}, {78, 43330}, {79, 48780}, {
          80, 53550}, {81, 60430}, {82, 66580}, {83, 75850}, {84, 87820}, {
          85, 94750}, {86, 102140}, {87, 109650}, {88, 118050}, {89, 123930}}

p1 = ListPlot[CostOfAdv, PlotStyle -> {PointSize[.02]}];

Linear Fit

c[x_] = Fit[CostOfAdv, {1, x}, x];

p0 = Plot[c[x], {x, 70, 89}, DisplayFunction -> Identity];

Show[p0, p1];

SSres = 
      (∑_(n = 1)^20)
          (((CostOfAdv[[n
                  ]][[2
                ]] - 
              c[CostOfAdv[[n
                    ]][[1
                  ]]]))^2)

meancost = 
      (1/20) (∑_(n = 1)^20)
             CostOfAdv[[n
              ]][[2
            ]]

SSreg = 
      (∑_(n = 1)^20)
          (((c[CostOfAdv[[n
                    ]][[1
                  ]]] - meancost))^2)

SStot = SSres + SSreg

(SSreg/SStot)

0.953070391452541780356`

Exponential Fit

Separate the cost data from the ordered pairs.

Cost = 
      Table[CostOfAdv[[n
            ]][[2
          ]], {n, 1, 20}]

We log the Cost data.

LoggedCost = N[Log[Cost]]

Recombine the ordered pairs.

LogCost = 
      Table[{CostOfAdv[[n
              ]][[1
            ]], 
          LoggedCost[[n]]}, {n, 1, 20}]

If this were true exponential growth, then this graph should be linear.

p2 = ListPlot[LogCost];

f[x_] = Fit[LogCost, {1, x}, x]

p3 = Plot[f[x], {x, 70, 89}, DisplayFunction -> Identity];

Show[p2, p3, DisplayFunction -> $DisplayFunction];

Transform back to original coordinates

g[x_] = (E^f)[x]

p4 = Plot[g[x], {x, 70, 89}, DisplayFunction -> Identity];

Show the model with the original data.

Show[p1, p4, DisplayFunction -> $DisplayFunction];

SSres = 
      (∑_(n = 1)^20)
          (((CostOfAdv[[n
                  ]][[2
                ]] - 
              g[CostOfAdv[[n
                    ]][[1
                  ]]]))^2)

meancost = 
      (1/20) (∑_(n = 1)^20)
             CostOfAdv[[n
              ]][[2
            ]]

SSreg = 
      (∑_(n = 1)^20)
          (((g[CostOfAdv[[n
                    ]][[1
                  ]]] - meancost))^2)

SStot = SSres + SSreg

(SSreg/SStot)

Power Fit

loggeddata = N[Log[CostOfAdv]]

p6 = ListPlot[loggeddata];

h[x_] = Fit[loggeddata, {1, x}, x]

p7 = Plot[h[x], {x, 70, 89}];

k[x_] = (E^h)[Log[x]]

p8 = Plot[k[x], {x, 70, 89}, DisplayFunction -> Identity];

Show[p1, p8];

Show[p0, p1, p4, p8, DisplayFunction -> $DisplayFunction];

SSres = 
      N[(∑_(n = 1)^20)
            (((CostOfAdv[[n
                    ]][[2
                  ]] - 
                k[CostOfAdv[[n
                      ]][[1
                    ]]]))^2)]

meancost = 
      (1/20) (∑_(n = 1)^20)
             CostOfAdv[[n
              ]][[2
            ]]

SSreg = 
      N[(∑_(n = 1)^20)
            (((k[CostOfAdv[[n
                      ]][[1
                    ]]] - meancost))^2)]

SStot = SSres + SSreg

(SSreg/SStot)

Fitting a Logistic Curve to Data:  Black Bear Data

This notebook fits a logistic curve to female black bear age/weight data.  The data is read from a graph in Smith and Clark, Black Bears in Arkansas:  Characteristics of a Successful Translocation, Journal of Mammology, Vol 75, No. 2, page 316.  The Mathematica analysis directly follows Fabio Cavallini's, Fitting a Logistic Curve to Data, College Mathematics Journal, Vol. 24, No. 3, May 1993.

The data where the first coordinate is the age class in years and the second coordinate is the weight in kilograms:

data={{1,35},{2,55},{3,68},{4,70},
{5,71},{6,75},{7,79},{8,82},{9,81},
{10,80},{11,78},{12,99},{13,99},{15,82}};

A scatter plot of the data:

ListPlot[data,PlotStyle->PointSize[.02],AxesLabel->{"Age Class(Years)","Weight(Kg)"},PlotRange->{20,100}];

Cavallini's article performs the analysis in a very elegant manner using vectors to deal with the data.  We follow the same method and notation.  See the text or Cavallini's article for details.  Splitting the data into two separate lists; one for age (T for time), one  for mass (Y).

T=Table[data[[n]][[1]],{n,1,14}]

Y=Table[data[[n]][[2]],{n,1,14}]

Create a vector of theoretical values in terms of parameters r and t0.  The result, H, will be a vector

H=1.0/(1.0+Exp[-r(T-t0)]);

The error vector is

error=(Y.Y)-(H.Y)^2/(H.H);

The key step in this analysis is chosing a wise starting point for the search algorithm.  Next we display a 3D plot of the error.    Since t0 must be within the data, choosing its range to be 0 to 16 is easy. A range for r may require some guesswork and repeated trials.

Plot3D[error, {r, 0.01, 1}, {t0, 0, 16}, 
      AxesLabel -> {"r",  "t0", "e    "}];

The next graph restricts the scale on the range to help identify the minimum.

Plot3D[error, {r, 0.01, 1}, {t0, 0, 16}, 
      AxesLabel -> {"r",  "t0", "e   "}, 
      PlotRange -> {0, 1000}];

We get a sense of the shape from this plot.  Next we will switch to a contour plot to flush out more details.  Recall that the contour plot is of error values for different r and t0 values.  The 3D Plot helps us determine the range of r and t0  to plot.  Darker regions indicate less error.

ContourPlot[error,{r,0.01,1},{t0,0,5}];

We will pick our initial guess from somewhere near the center of the black region.  The following gives the error for these parameter choices.

error/.{r->.4,t0->1}

Use this guess to start the FindMinimum algorithm.  You can experiment with other r, t0 combinations to try to get a smaller error, but all we need is close enough for the search routine to work correctly.

FindMinimum[error,{r,0.4},{t0,1}]

The minimum error is 585.75 at parameter values r=.441 and t0=1.152

Find K for these parameter values, using Cavallini's method.

K=((H.Y)/(H.H))/.{r->0.440661014003951072466`,t0->1.15244778447574305923`}

Another copy of the scatterplot of the data.

FigData=ListPlot[Transpose[{T,Y}],PlotStyle->PointSize[0.02],PlotRange->{20,100}];

A plot of the fitted (found parameters inserted) logistic curve.

FigModel=Plot[K/(1+Exp[-0.440661014003951072466`(t-1.15244778447574305923`)]),
{t,Min[T],Max[T]},PlotRange->{20,100}];

Finally the theoretical curve displayed with the data.

Show[FigData,FigModel,AxesLabel->{"Age Class(Years)","Weight(Kg)"},DisplayFunction->$DisplayFunction];

Fitting the Lynx Data

Data from http://biology.wsc.mass.edu/biology/courses/hoag/popmodel/fall99/homework/empirical/lynx.htm, who took it from our textbook!

\(lynxdata
       = {{1919.0, 8378.0}, 
{1920.0, 6456.0}, 
{1921.0, 11617.0}, 
        
{1922.0, 17202.0}, 
{1923.0, 26381.0}, 
{1924.0, 29529.0}, 
        
{1925.0, 33027.0}, 
{1926.0, 28619.0}, 
{1927.0, 21363.0}, 
        
{1928.0, 11582.0}, 
{1929.0, 7572.0}, 
{1930.0, 7957.0}, 
        
{1931.0, 8410.0}, 
{1932.0, 11916.0}, 
{1933.0, 16781.0}, 
        
{1934.0, 22012.0}, 
{1935.0, 22448.0}, 
{1936.0, 17534.0}, 
        
{1937.0, 10523.0}, 
{1938.0, 8079.0}, 
{1939.0, 7411.0}, 
        
{1940.0, 6583.0}, 
{1941.0, 6979.0}, 
{1942.0, 7544.0}, 
        
{1943.0, 10164.0}, 
{1944.0, 12259.0}, 
{1945.0, 9306.0}, 
        
{1946.0, 8129.0}, 
{1947.0, 6548.0}, 
{1948.0, 4083.0}, 
        
{1949.0, 3714.0}, 
{1950.0, 9592.0}, 
{1951.0, 6653.0}, 
        
{1952.0, 12636.0}, 
{1953.0, 10876.0}, 
{1954.0, 13876.0}, 
        
{1955.0, 9660.0}, 
{1956.0, 8397.0}, 
{1957.0, 8958.0}
}; 


    lynxPlot\  = \ 
      ListPlot[lynxdata, AxesLabel -> {"Year", "Return"}, 
        PlotStyle -> {PointSize[.02]}]; \)

linearlynx[x_] = Fit[lynxdata, {1, x}, x]nonlinearlynx[x_] = 
      Fit[lynxdata, {1, x, Sin[2  Pi x/10], Cos[2  Pi x/10]}, x]

\(linearlynxext\  = \ 
      Plot[linearlynx[x], {x, 1919, 1957}, 
		
        AxesLabel -> {"Year", "Return"}, 
        DisplayFunction -> Identity]; 

    nonlinearlynxext\  = \ 
      Plot[nonlinearlynx[x], {x, 1919, 1957}, 
        AxesLabel -> {"Year", "Return"}, 
        DisplayFunction -> Identity]; 

    Show[lynxPlot, linearlynxext, \ nonlinearlynxext, 
      DisplayFunction -> $DisplayFunction]; \)

detrendedlynx = 
      Table[{lynxdata[[n
              ]][[1
            ]], 
          lynxdata[[n
                ]][[2
              ]] - 
            linearlynx[
              lynxdata[[n
                  ]][[1
                ]]]},  {n, 1, Length[lynxdata]}]detrendedlynxPlot  =  
      ListPlot[detrendedlynx, AxesLabel -> {"Year", "Return"}, 
        PlotStyle -> {PointSize[.02]}];

I'm going to do the detrending a little differently than they do: I'll start with an exponentially decaying function to model the general decline, rather than a linear function. To do this, we log-transform the data, then fit the linear function to the transformed data; this is followed by a back-transform:

loglynx2 = 
      Table[{lynxdata[[n
              ]][[1
            ]], 
          Log[lynxdata[[n
                ]][[2
              ]]]},  {n, 1, Length[lynxdata]}]expolynx[x_] = Fit[loglynx2, {1, x}, x]trend[x_] = (E^expolynx)[x]

We then detrend the data, using the exponential trend function:

detrendedlynx2 = 
      Table[{lynxdata[[n
              ]][[1
            ]], 
          lynxdata[[n
                ]][[2
              ]] - 
            trend[lynxdata[[n
                  ]][[1
                ]]]},  {n, 1, Length[lynxdata]}]detrendedlynx2Plot  =  
      ListPlot[detrendedlynx2, AxesLabel -> {"Year", "Return"}, 
        PlotStyle -> {PointSize[.02]}];

Now we compute the envelope in the same way that they did it in the book:

envelope = 
      Table[{detrendedlynx2[[5  n + 2]][[1]], 
          Abs[detrendedlynx2[[5  n + 2]][[2]]]}, {n, 0, 7}]envelopePlot  =  
      ListPlot[envelope, AxesLabel -> {"Year", "Return"}, 
        PlotStyle -> {PointSize[.02]}];

We use the same procedure, of computing an exponential model by log-transforming the envelope data, and then backtransforming the linear model:

logenvelope = 
      Table[{envelope[[n
              ]][[1
            ]], 
          Log[envelope[[n
                ]][[2
              ]]]},  {n, 1, Length[envelope]}]expoenvelope[x_] = Fit[logenvelope, {1, x}, x]trend2[x_] = (E^expoenvelope)[x]

We compute the non-linear model to the by using only the detrended data (plus we introduce a constant term).

nonlineardetrended2[x_] = 
      Fit[detrendedlynx2, {1, trend2[x]*Sin[2  Pi x/10], 
          trend2[x]*Cos[2  Pi x/10]}, x]

Finally, we compute the general model by combining the two models: the general exponential decay trend, and the detrended data model:

model[x_] = nonlineardetrended2[x] + trend[x]\(modelext\  = \ 
      Plot[model[x], {x, 1919, 1957}, 
		
        AxesLabel -> {"Year", "Return"}, 
        DisplayFunction -> Identity]; 

    Show[lynxPlot, modelext, DisplayFunction -> $DisplayFunction]; \)

Exercise 3

Data from Exercise #3, p. 226

trial
       = {{0, 2.2},  {1,  5.2},  {2,  4.8},  {3,  11.1},  {4, 13.9}, 
         {5, 17.1}}\(Needs["Statistics`DescriptiveStatistics`"]; 

    meanx = \(∑\+\(n = 1\)\%\(Length[trial]\)\(trial[[n
              ]]\)[[1
            ]]\)\/Length[trial]\)meany = 
      (((∑_(n = 1)^(Length[trial]))trial[[n
              ]][[2
            ]])/Length)[trial]\(xs = Table[\(trial[\([n]\)]\)[\([1]\)], {n, 1, 6}]; 

    ys = Table[\(trial[\([n]\)]\)[\([2]\)], {n, 1, 6}]; \)

cov = 
      1/(Length[trial] - 1)*
        (∑_(n = 1)^(Length[trial]))(xs[[n]] - meanx)*
            (ys[[n]] - meany)beta1 = cov/Variance[xs]*1.0

StandardDeviation[xs]*1.0StandardDeviation[ys]r = cov/StandardDeviation[xs]/StandardDeviation[ys]*1.0

Ex3[x_] = Fit[trial, {1, x}, x]

Next we fit a line to the first season data and overlay the line and the data.

trialPlot = ListPlot[trial, AxesLabel->{"x","y"},PlotStyle->{PointSize[.02]},DisplayFunction->Identity];

trialext = Plot[Ex3[x],{x,0,5},AxesLabel->{"x","y"},DisplayFunction->Identity];

Show[trialPlot,trialext,DisplayFunction->$DisplayFunction];

data = trial;

SSres = 
      (∑_(n = 1)^(Length[data]))(((
            data[[n
                  ]][[2
                ]] - 
              Ex3[data[[n
                    ]][[1
                  ]]]))^2)

mean = 
      (((∑_(n = 1)^(Length[data]))data[[n
              ]][[2
            ]])/Length)[data]

SSreg = 
      (∑_(n = 1)^(Length[data]))(((
            Ex3[data[[n
                    ]][[1
                  ]]] - mean))^2)

SStot = SSres + SSreg

(SSreg/SStot)