<< Graphics`

Curve Fitting and Population Growth

This lab takes you through a very simple-minded approach to fitting different curves to some population data and examining the results.  This is in no way meant to be a complete statistical treatment of this subject.  Mathematica offers much more sophisticated tools to fit curves and analyze the results than what I use in this lab.  I chose these because they do not require any specialized knowledge of statistics to use.

The data

For this lab, I went to the US Census Bureau website and downloaded population summary data for the United States from 1950-1980.  I saved it in a "comma-delimited" text file (a csv file) and read it into Mathematica using the following command (you will get an error message if you execute it):

populationData = 
      Import["c:\TEMP\population_data_thru_1980.csv"]

So you don't have to worry about importing the data, you should execute the following statement to define our population data:

populationData = {{1950, 152271000}, {1951, 154878000}, {1952, 
          157553000}, {1953, 160184000}, {1954, 163026000}, {1955, 
          165931000}, {1956, 168903000}, {1957, 171984000}, {1958, 
          174882000}, {1959, 177830000}, {1960, 180671000}, {1961, 
          183691000}, {1962, 186538000}, {1963, 189242000}, {1964, 
          191889000}, {1965, 194303000}, {1966, 196560000}, {1967, 
          198712000}, {1968, 200706000}, {1969, 202677000}, {1970, 
          205052000}, {1971, 207661000}, {1972, 209896000}, {1973, 
          211909000}, {1974, 213854000}, {1975, 215973000}, {1976, 
          218035000}, {1977, 220239000}, {1978, 222585000}, {1979, 
          225055000}, {1980, 227726463}}

(I only included the Import function in case some of you want to experiment with importing your own data.)

Let's take a look at what this data looks like:

populationPlot = 
      ListPlot[populationData, PlotRange -> {0, Max[populationData]}, 
        PlotStyle -> {Red, PointSize[0.01]}]

Notice that we use the ListPlot function to plot discrete data like this.  I explicitly specify the PlotRange to be sure I see everything and I use PlotStyle to set the color and size of the dots.

Linear Fit

This looks basically linear, so let's define our linear least-square fit using the Mathematica Fit function.  We tell it to fit this data to the functions 1 and t (a constant and t to the first power):

linearFit[t_] = Fit[populationData, {1, t}, t]

Let's graph this and see what it looks like (I use the same PlotRange as before, but make it blue):

linearPlot = 
      Plot[linearFit[t], {t, 1950, 1980}, 
        PlotRange -> {0, Max[populationData]}, 
        PlotStyle -> Blue]

And now, display them together:

Show[linearPlot, populationPlot]

This doesn't look like too bad a fit.

Quadratic Fit

On the other hand,perhaps we could do better with a quadratic (parabolic) fit:

quadraticFit[t_] = Fit[populationData, {1, t, (t^2)}, t]

quadraticPlot = 
      Plot[quadraticFit[t], {t, 1950, 1980}, 
        PlotRange -> {0, Max[populationData]}, 
        PlotStyle -> Blue]

Show[quadraticPlot, populationPlot]

This one looks even better.

Exponential Fit

Of course, in almost all of your math courses, you have probably heard that population growth can best be modeled by an exponential growth function.  This makes sense (or we have told you it does), but is it true?

Unfortunately, fitting data to an exponential growth equation y = (Subscript[y, 0]) (e^(k t)) takes a little more work in Mathematica.  To do this, first notice that \(TraditionalForm\`y = \(y\_0\ e\^\(k\ t\) = \(\(e\^\(\(\ \)\(a\)\)\) e\^\(\(\ \)\(k\ t\)\) = e\^\(\(\ \)\(a + k\ t\)\)\)\)\) for some constant a = ln (Subscript[y, 0]).  Thus, \(TraditionalForm\`ln\ y = \(ln\ e\^\(\(\ \)\(a + k\ t\)\) = a + k\ t\)\).  So, in other words, the natural log of y is a line for exponential growth.  

So, our technique will be the following:
1)  Take the natural log of all the y values (the population).
2)  Fit those log's to a line, using the Fit function.
3)  Exponentiate that line to get back to population (i.e., to "cancel out" the log's we took in part 1).

Take the log of all the population data

1) Taking the natural log of all of our population data is a little tricky.We want to define a new data set whose first coordinate for each point is the year,but whose second coordinate itLog(population)(a reminder:Mathematicauses Log for the natural log,ln).In a list of data,you can access different elements of that list using the "double-bracket" notation:

populationData[[1]]

populationData[[2]]

What we really need to access is the set of all populations.  One way to do this is to take the Transpose of the populationData:

Transpose[populationData]

Now, the first element contains all the years:

Transpose[populationData][[1]]

And the second element contains all the populations (in the correct order):

Transpose[populationData][[2]]

One particularly cool trick in Mathematica is that you can take the natural log of all of these at one time:

Log[Transpose[populationData][[2]]]

Of course, after we finish, we need to Transpose the data back into the original form (so each element contains a year and a population).  This can all be combined together into the following command:

logPopulationData = 
      Transpose[{Transpose[populationData][[1]], 
            Log[Transpose[populationData][[2]]]}] // N

Fit the log's to a line

2)  Fit the log's of the population to a line.  This is easy:

logLine = Fit[logPopulationData, {1, t}, t]

Eponentiate the result

3)  And finally, we exponentiate the result to get back to the actual population (instead of the log of the population):

exponentialFit[t_] = (E^logLine)

exponentialPlot = 
      Plot[exponentialFit[t], {t, 1950, 1980}, 
        PlotRange -> {0, Max[populationData]}, 
        PlotStyle -> Yellow]

Show[exponentialPlot, populationPlot]

This doesn't seem terribly unreasonable either.

Compare the different fits

So, which one of these three methods actually gives the best results?  Let's view them all together on one graph:

Show[linearPlot, quadraticPlot, exponentialPlot, 
      populationPlot]

Compare the projections to actual data

In case you were wondering why I only took data from 1950-1980, it is because I want to now load data for 1950-2000 and see which of our projections (if any) predicts the actual population growth accurately.  Let's load the complete data now:

populationData2 = {{1950, 152271000}, {1951, 154878000}, {1952, 
          157553000}, {1953, 160184000}, {1954, 163026000}, {1955, 
          165931000}, {1956, 168903000}, {1957, 171984000}, {1958, 
          174882000}, {1959, 177830000}, {1960, 180671000}, {1961, 
          183691000}, {1962, 186538000}, {1963, 189242000}, {1964, 
          191889000}, {1965, 194303000}, {1966, 196560000}, {1967, 
          198712000}, {1968, 200706000}, {1969, 202677000}, {1970, 
          205052000}, {1971, 207661000}, {1972, 209896000}, {1973, 
          211909000}, {1974, 213854000}, {1975, 215973000}, {1976, 
          218035000}, {1977, 220239000}, {1978, 222585000}, {1979, 
          225055000}, {1980, 227726463}, {1981, 229966237}, {1982, 
          232187835}, {1983, 234307207}, {1984, 236348292}, {1985, 
          238466283}, {1986, 240650755}, {1987, 242803533}, {1988, 
          245021414}, {1989, 247341697}, {1990, 250131894}, {1991, 
          253492503}, {1992, 256894189}, {1993, 260255352}, {1994, 
          263435673}, {1995, 266557091}, {1996, 269667391}, {1997, 
          272911760}, {1998, 276115288}, {1999, 279294713}, {2000, 
          282338631}}

populationPlot2 = 
      ListPlot[populationData2, PlotRange -> {0, Max[populationData2]}, 
        PlotStyle -> {Red, PointSize[0.01]}]

Let's plot this new data along with our three different functions (linearFit, quadraticFit, and exponentialFit) together now.  Notice that I need to extend the domains on the functions to graph this (keeping the same colors I used above):

Show[
      Plot[{linearFit[t], quadraticFit[t], exponentialFit[t]}, {t, 1950, 
          2000}, PlotRange -> {0, Max[populationData2]}, 
        PlotStyle -> {Green, Blue, Yellow}], populationPlot2]

Notice that these three functions now begin to diverge from each other significantly.

Questions

Which of our functions looks to be the "best"?  Why?  Do you expect this to be true as even more data is collected for future years?

Graph the three functions through the year 2050 and discuss what sort of difference it makes which of our functions you choose to use to model the population growth.

A longer baseline:  Notice that our data started in 1950.  What happens if we use data for a longer period of time?  Use the data below which includes data for every 10 years between 1900 and 2000 and compare it with our three functions.  Which is best?  Why?

populationData3 = {{2000, 281421906}, {1990, 248709873}, {1980, 
          226542199}, {1970, 203302031}, {1960, 179323175}, {1950, 
          151325798}, {1940, 132164569}, {1930, 123202624}, {1920, 
          106021537}, {1910, 92228496}, {1900, 76212168}}

Use the new data (populationData3) to fit new linear, quadratic, and exponential functions.  Discuss which of these three functions works best (and why).

Compare the new functions to the more detailed data in populationData2 and discuss whether the new functions fit that data as well.

So, can the population really be modeled by an exponential function or not?  What other improtant observations can you make about curve-fitting and data from this exercise?