<Text-field style="Heading 1" layout="Heading 1">Modeling Change: Springs, Driving Safety, Radioactivity, Trees, Fish, and Mammals Note: You may notice differences between this Maple worksheet and the equivalent Mathematica notebook. These differences were introduced to preserve the content of these modules and were necessary because of major functional differences between Maple and Mathematica.</Text-field>

<Text-field style="_cstyle307" layout="Heading 1">Introduction</Text-field>

OBJECTIVE: Practice a modeling process: consider a behavior, observe data, fit a model, analyze the error, improve the model if appropriate, interpret the model, and make predictions.

Mathematical models help us better understand things in the world around us, like the behavior of springs, safe driving practices, radioactivity in medicine, the growth of trees and fish, and the biology of mammals. In this module you will learn how to construct mathematical models, how analyze and improve them, and how to use them to study the behavior you are modeling, to make predictions. As you will see, a computer algebra system like Maple is a powerful tool that will aid you in your mathematical modeling. In this worksheet, we use the two Maple packages plots and stats. We initialize them at the beginning of Part I.

<Text-field style="_cstyle298" layout="Heading 1">Technology Guidelines</Text-field>

NOTE: If you have just finished a worksheet, restart Maple before executing a new worksheet. TO OPEN SECTIONS, Click on the PLUS sign at the left hand side of the screen or select Expand All Sections from the View drop down menu.

TO STOP AN EXECUTION Click on STOP button from the toolbar.

ORDER OF EXECUTION Execute commands in the order given. Do not skip any Maple Input lines within a given worksheet

Alternatively, you can execute the entire worksheet by selecting the Execute Worksheet command from the Edit drop down menu.

SAVING WORKSHEETS.

You can save anytime to any directory you choose, and it is wise to save often. EXPERIENCING MAJOR PROBLEMS Save if appropriate, and then shut down Maple and start it up again.

<Text-field style="Heading 1" layout="Heading 1">Part I: Spring Elongation</Text-field>

<Text-field style="_cstyle257" layout="Bullet Item">Observing the Data</Text-field>

The response of a spring to various loads can be modeled in order to design a vehicle such as a tank, utility vehicle, or luxury car that responds to road conditions in a desired way. (See also "Bungee Cord Jumping: A Classroom Experiment," another module in this supplement.) We conducted an experiment to measure the stretch of a spring in inches as a function of the number of units of mass placed on the spring. The following list includes these data with the first element in each ordered pair being the mass on the spring and the second element its corresponding stretch.

restart; with(plots): with(stats):

data:= [[0, 0], [1, .875], [2, 1.721], [3, 2.641], [4, 3.531], [5, 4.391],[6, 5.241], [7, 6.120], [8, 6.992], [9, 7.869], [10, 8.741]];

matrix([[`mass `,`stretch (in)`],op(data)]);

Next, we plot the data to see if there is a recognizable pattern and name the plot p1 for later use.

xlabel:="mass": ylabel:="stretch (in)": p1:=pointplot(data, labels=[xlabel,ylabel], labeldirections=[HORIZONTAL, VERTICAL]): print(p1);

<Text-field style="_cstyle265" layout="Bullet Item">Designing a Model</Text-field>

The data strongly suggest a linear relationship. Now we use the fit[leastsquare]( ) function to find the line of best-fit or the linear regression line.

data2 := [[seq(data[i,1],i=1..nops(data))],[seq(data[i,2],i=1..nops(data))]]: eqn:=fit[leastsquare[[x,y],y=a*x]]([data2[1],data2[2]]);

The constant of proportionality is 0.8747. Next, we plot the regression line, save it as p2, and then show the plots of the data and the regression line together on the same graph.

p2:=implicitplot(eqn , x=0..10, y=-10..10, labels=[`mass`, `stretch (in)`], labeldirections=[HORIZONTAL, VERTICAL], numpoints=100): print(p2); print(display({p1,p2}));

<Text-field style="_cstyle266" layout="Bullet Item">Assessing the Errors</Text-field>

As we expected, the line appears to fit the data very well. We can provide further verification by calculating the residual errors for each data point and plotting them. First, we need to use the model to calculate the values of stretch for each mass value in the original data set.

f:=unapply(rhs(eqn),x): predictvalues:=[seq([i, f(i)], i=0..10)];

The residual error (or residual) is the difference between the measured value and the value the model predicts for each amount of mass.

res:=data-predictvalues: residuals:=[seq([i-1, res[i,2]],i=1..11)];

pointplot(residuals, labels=["mass","residual error"], labeldirections=[HORIZONTAL,VERTICAL],view=[-2..12,-0.03..0.03]);

The residuals show that the differences are indeed small. Small is a relative term, and sometimes it is more meaningful to look at the error in relation to the size of the quantity being estimated. For all but the first data point, we calculate the relative error by dividing the measured value into the residual and multiplying by 100 to express this relative error as a percentage. Then, we plot the percentage error for each data point.

percenterrors:=[seq([i-1, 100*residuals[i,2]/data[i,2]],i=2..11)];

pointplot(percenterrors, labels=["mass","% error"], labeldirections=[HORIZONTAL,VERTICAL],view=[-2..12,-2..1]);

The maximum relative error for the linear model is LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JS1JI21vR0YkNi1RKiZ1bWludXMwO0YnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGOi8lKXN0cmV0Y2h5R0Y6LyUqc3ltbWV0cmljR0Y6LyUobGFyZ2VvcEdGOi8lLm1vdmFibGVsaW1pdHNHRjovJSdhY2NlbnRHRjovJSdsc3BhY2VHUSwwLjIyMjIyMjJlbUYnLyUncnNwYWNlR0ZJLUkjbW5HRiQ2JFEkMS43RidGNUY1RitGNQ==%. Not bad!

<Text-field style="Heading 1" layout="Heading 1">Part II: Safe Following Distance</Text-field>

<Text-field style="_cstyle267" layout="Bullet Item">Observing the Data</Text-field>

A rule of thumb that is often given for safe following distance is to allow 2 seconds between your car and the car in front of you. Is this rule reasonable? The following data set contains ordered pairs of values. The first element of each ordered pair is the traveling speed, and the second element is the total stopping distance, the distance traveled by the car during the driver's reaction time plus the distance required for the vehicle to come to a stop with full braking.

data := [[20., 22 + 32], [25, 28 + 47], [30, 33 + 65], [35, 39 + 87], [40, 44 + 112], [45, 50 + 140], [50, 55 + 171], [55, 61 + 204], [60, 66 + 241], [65, 72 + 282], [70, 77 + 325], [75, 83 + 376]];

matrix([[`speed`, `stopping dist`],op(data)]);

We plot the data to see if there is a recognizable pattern

xlabel:="speed (mph)": ylabel:="stopping distance (ft)": p1:=pointplot(data, labels=[xlabel,ylabel], view=[-5..75,-50..500],labeldirections=[HORIZONTAL, VERTICAL]): print(p1);

<Text-field style="_cstyle268" layout="Bullet Item">Designing a Model</Text-field>

The data suggest a possible quadratic relationship. Now we use the fit[leastsquare]( ) function to find the curve of best-fit or the regression curve.

temp := [[seq(data[i,1],i=1..nops(data))],[seq(data[i,2],i=1..nops(data))]]: eqn:=fit[leastsquare[[x,y], y=a+b*x+c*x^2] ](temp); f:=unapply(rhs(eqn),x):

Next, we plot the regression curve and then show the plot of the data and the model together on the same graph.

p2:=plot(f(x), x=0..85, labels=[xlabel,ylabel],view=[-5..75,-50..500], labeldirections=[HORIZONTAL, VERTICAL]): print(p2);

print(display({p1,p2}));

The best-fit quadratic function appears to fit the data very well.

<Text-field style="_cstyle269" layout="Bullet Item">Assessing the Errors</Text-field>

We can provide further verification by calculating the residual errors for each data point and plotting them. First, we calculate the stopping distances predicted by the model, and then we calculate the residual errors and plot them.

f:=unapply(rhs(eqn),x): predictvalues:=[seq([i*5., f(i*5.)], i=4..15)];

matrix([[`speed `,`predicted values`],op(predictvalues)]);

res:=data-predictvalues: residuals:=[seq([20+5*(i-1), res[i,2]],i=1..nops(data))];

pointplot(residuals, labels=[xlabel, "residual error"], labeldirections=[HORIZONTAL, VERTICAL]);

Once again we should look at the error in relation to the size of the quantity being estimated. For all but the first data point, we can calculate the relative error as follows.

percenterrors:=[seq([20+5*(i-1), 100*residuals[i,2]/data[i,2]],i=2..11)];

pointplot(percenterrors, labels=[xlabel, "percentage error"], labeldirections=[HORIZONTAL, VERTICAL],view=[-5..70,-0.5..0.6]);

The maximum relatvie error of this model is about 0.6%.

<Text-field style="_cstyle270" layout="Bullet Item">Improving the Model</Text-field>

As indicated at the beginning of this module, a rule of thumb that is often given for safe following distance is to allow 2 seconds between your car and the car in front of you. To use this rule, you would note when the car in front of you passes some marker on or near the roadway, and then you count the time it takes you to get to the marker. This time should be at least 2 seconds. The rationale behind this rule is that if the car in front of you were suddently to come to a complete stop, the 2-second separation distance would give you enough time to stop, to avoid hitting the vehicle in front of you. To test the 2-second rule of thumb, we will calculate how far your car travels in 2 seconds at various speeds and then compare these distances with the corresponding stopping distances in our data set.

The first entry in each element of the following table is the traveling speed in miles per hour, the second entry is the 2-second distance between the two cars in feet, and the third entry is the stopping distance in feet, taken from the test data. (To calculate the distance traveled in 2 seconds, we convert the traveling speeds from mph to ft/sec by using the conversion, 60 mph = 88 ft/sec.)

data2 := [seq([data[i, 1], (data[i, 1]*88./60.0)*2, data[i, 2]], i=1..12)];

matrix([[`v (mph)`,`2-sec dist. (ft)`, `stopping dist(ft)`],op(data2)]);

The data in the table show that for speeds greater than 20 mph, the stopping distance exceeds the 2-second separation distance.

Let's find an improved model to provide the basis for a better rule of thumb. To do this, we will look at the problem the other way around. We know the stopping distance for various speeds. If we force the stopping distance and the separation distance to be equal, then we can calculate separation time for each speed. We do this by taking the stopping distance (in feet) and dividing it by the travel speed (in ft/sec). In the following table, the first entry in each element is the traveling speed in miles per hour, and the second entry is the separation time in seconds that is required to ensure a separation distance equal to the stopping distance.

data3:=[seq([data[i, 1], data[i, 2]/( data[i, 1]*88./60)], i=1..(nops(data)))];

matrix([[`speed (mph)`,`separation time (sec)`],op(data3)]);

Next time plot the separation time versus speed.

xlabel:=`speed (mph)`: ylabel:=`separation time (sec)`: p1:=pointplot(data3, labels=[xlabel,ylabel], labeldirections=[HORIZONTAL, VERTICAL], view=[-5..80,-1..5]): print(p1);

It appears that the separation time is a linear function of the travel speed and not a constant function as the 2-second rule suggests. Let's find the best-fit linear function for the separation-time versus speed data and plot it together with the data points.

temp:= [[seq(data3[i,1],i=1..nops(data3))],[seq(data3[i,2],i=1..nops(data3))]]: fit[leastsquare[[x,y], y=a+b*x] ](temp): f:=unapply(rhs(%),x);

p2:=plot(f(x), x=0..85, labels=[xlabel, ylabel], labeldirections=[HORIZONTAL, VERTICAL], view=[-5..80,-1..5]): print(p2);

print(display({p1,p2}));

<Text-field style="Heading 1" layout="Heading 1">You Try It: Make a Better Rule of Thumb</Text-field>

Using the Model

Based upon the results of the analysis in Part II, formulate a better rule of thumb for separation time versus travel speed. Keep in mind that a rule of thumb should be easy to remember, easy to use, and most importantly, it should be correct.

<Text-field style="Heading 1" layout="Heading 1">Part III: Radioactivity</Text-field>

<Text-field style="_cstyle272" layout="Bullet Item">Observing the Data</Text-field>

A radioactive dye is injected into a patient's veins to facilitate an X-ray procedure. Measuring the radioactivity in counts per minute every minute for 10 minutes yielded the table of values shown below.

data := [[0, 10023], [1, 8174], [2, 6693], [3, 5500], [4, 4489], [5, 3683], [6,3061], [7, 2479], [8, 2045], [9, 1645], [10, 1326]];

matrix([[`Time (min)`,`Radioactivity (cpm)`],op(data)]);

We begin by plotting the data points.

xlabel := "x (min)": ylabel := "y (cpm)": p1:=pointplot(data, labels=[xlabel, ylabel], labeldirections=[HORIZONTAL, VERTICAL], view=[-2..10,-1000..11000]): print(p1);

<Text-field style="_cstyle273" layout="Bullet Item">Designing a Model</Text-field>

There appears to be a trend that we can capture with a mathematical model, and now we try to find a suitable model using the fit[leastsquares]( ) function. What does the function look like to you? A decaying exponential? We try to find a function of the form LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkjbWlHRiQ2I1EhRictRiM2KC1GLDYmUSJ5RicvJSVib2xkR1EmZmFsc2VGJy8lJ2l0YWxpY0dRJXRydWVGJy8lLG1hdGh2YXJpYW50R1EnaXRhbGljRictSSNtb0dGJDYuUSI9RidGNC9GO1Enbm9ybWFsRicvJSZmZW5jZUdGNi8lKnNlcGFyYXRvckdGNi8lKXN0cmV0Y2h5R0Y2LyUqc3ltbWV0cmljR0Y2LyUobGFyZ2VvcEdGNi8lLm1vdmFibGVsaW1pdHNHRjYvJSdhY2NlbnRHRjYvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZTLUYjNictRiw2JlEiYUYnRjRGN0Y6LUY+Ni1RMSZJbnZpc2libGVUaW1lcztGJ0ZBRkNGRUZHRklGS0ZNRk8vRlJRJjAuMGVtRicvRlVGaW4tSSVtc3VwR0YkNiUtRj42LlEvJkV4cG9uZW50aWFsRTtGJ0Y0RkFGQ0ZFRkdGSUZLRk1GT0Zobi9GVVEsMC4xMTExMTExZW1GJy1GIzYnLUY+Ni5RKiZ1bWludXMwO0YnRjRGQUZDRkVGR0ZJRktGTUZPL0ZSUSwwLjIyMjIyMjJlbUYnL0ZVRmlvLUYjNictRiw2JlEia0YnRjRGN0Y6RmVuLUYsNiZRInhGJ0Y0RjdGOi9GO1ElYm9sZEYnLyUrZm9udHdlaWdodEdGZHBGK0ZjcEZlcC8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGY3BGZXBGK0ZjcEZlcEYrRmNwRmVw, where we vary the value of k and the computer finds the a to find the best-fit function of this form.

temp:= [[seq(data[i,1],i=1..nops(data))],[seq(data[i,2],i=1..nops(data))]]: fit[leastsquare[[x,y], y=a*exp(-.1*x),{a}]](temp): f2:=unapply(rhs(%),x);

Next, we plot our model function, save it as p2, and then show the plots of the data and the model function together on the same graph.

p2:=plot(f2(x), x=0..12, y=0..10000, labels=["x (min)","y (cpm)"], labeldirections=[HORIZONTAL, VERTICAL], view=[-2..10,-1000..11000]): print(p2);

print(display({p1,p2}));

<Text-field style="_cstyle274" layout="Bullet Item">Assessing the Errors</Text-field>

Our model doesn't appear to be very good! Let's quantify just how good (or bad) it is. You guessed right we should calculate the residual errors to do this quantification.

predictvalues:=[seq([i, f2(i)], i=0..10)]: res:=data-predictvalues: residuals:=[seq([i-1, res[i,2]],i=1..11)];

To get an overall measure of the error for all of the data points in the set, you might be inclined to calculate the average (mean) all of the individual or local residuals, but this can be very misleading.

residavg:=add(residuals[i,2],i=1..11)/nops(residuals);

As you can see, the mean residual error is small relative to the average size of the measured values of radioactivity in the data set. On this basis, we might be led to believe that our model isn't so bad after all. Wrong! The problem is that the individual or local errors are actually sizable when compared to the measured radioactivity counts, but because some of them are positive and some are negative, they tend to cancel each other when we add them together to calculate their mean value. Look at the values in the list of residual errors or look at the vertical differences between the data points and the fit function on the graph, and it is easy to see that the average of the residual errors is misleading. There is an infinite number of ways to address this canceling problem, but one of the most common is to calculate the sum of the squares of the residuals and try to find the smallest or least value of this sum of squared residuals, hence we have the term "least squares." Squaring the residuals removes the canceling effect that occurs when we add them together for a measure of the global error. The Maple fit[leastsquare]( ) function uses least squares, and some variations of it, to find a best-fit function for a set of data. Finding the minimum value for the sum of squared residuals is a problem that can be solved using calculus, and you will study this problem later on, but for now you can do it by trial and error. Let's get back to our radioactivity problem. What we do now is to calculate the sum of the squares of the residuals for our data set.

residsquaresum:=add((residuals[i,2])^2,i=1..11);

Before we compare this value, we calculate the mean of the squared residuals.

msresiduals:=residsquaresum/nops(residuals);

Comparing the average of the squared residuals with the average of the radioactivity counts in the data set would be like comparing apples with oranges because msresiduals is an average of squares, whereas the average radioactivity count is not an average of squared values. To make a fair comparison, we take the square root of the mean of the squares.

rmsresiduals:=sqrt(msresiduals);

The value that we calculate in the preceding step is the root of the mean of the squares of the residuals and is oftentimes called the root-mean-square or "rms" value of the residuals.

<Text-field style="Heading 1" layout="Heading 1">You Try It: Improving the Model and Making Predictions</Text-field>

Improving the Model

In Part III, we calculated the "rms" value of the residual errors, a number that we can use for a fair comparison with the average of the measured radioactivity values. As we initially expected, this comparison leads us to conclude that our model needs improvement. We will leave that up to you, but, to help out, we group all of the commands that you need for the error analysis into one cell, the one that follows. Use trial and error to find a better model function by changing the value of k to reduce (possibly minimize) the sum of the squared residuals (which will also minimize the "rms" value).

A comment about errors in modeling: Keep in mind that in mathematical modeling you may not be able to completely eliminate errors. This is because of random errors that occur in measurements due to the limited precision of all measuring instruments, and because of systematic errors (i.e., errors that follow a pattern) due to shortcomings of the model and/or possible defects in the measuring tools. Systematic errors can often be reduced or eliminated by refining the model and repairing or recalibrating the measuring equipment, but random errors are still unavoidable. We can reduce the magnitudes of random errors by using more precise instruments, but we cannot eliminate them. The errors may be smaller, but they are always present.

k:= 0.2: temp:= [[seq(data[i,1],i=1..nops(data))],[seq(data[i,2],i=1..nops(data))]]: fit[leastsquare[[x,y], y=a*exp(-k*x),{a}]](temp); f:=unapply(rhs(%),x): p2:=plot(f(x), x=0..10, labels=[xlabel, ylabel]): print(display({p2, p1})); residuals := [seq([data[i, 1], data[i, 2] - f(data[i, 1])], i=1..nops(data))]: residsquaresum:=add((residuals[i,2])^2,i=1..11): mresiduals:=residsquaresum/nops(data): rmsresiduals:=sqrt(mresiduals);

Making Predictions

Now use your improved model to determine when the radioactivity will fall below 500 counts per minute.

solve(f(x)=500,x);

Your improved model should predict that the radioactivity will be below 500 cpm after about 15 minutes. Is that what you get?

<Text-field style="Heading 1" layout="Heading 1">Part IV: Ponderosa Pines</Text-field>

<Text-field style="_cstyle276" layout="Bullet Item">Observing the Data</Text-field>

In the table that follows, the girth of the pine tree (the distance arond the tree at shoulder height) is measured in inches, and the volume of usable lumber obtained from the tree is measured in board feet (bf). We will formulate and test the following models: that usable board feet is proportional to (a) the square of the girth and (b) the cube of the girth. Which is better? Does one model provide a better explanation than the other?

data := [[17, 19], [19, 25], [20, 32], [23, 57], [25, 71], [28, 113], [32,123], [38, 252], [39, 259], [41, 294]];

matrix([[`girth (in)`,`lumber (bf)`],op(data)]);

Next, we plot the data to see if there is a recognizable pattern.

xlabel:=`girth(in)`: ylabel:=`lumber (bf)`: p1:=pointplot(data, labels=[xlabel,ylabel], labeldirections=[HORIZONTAL, VERTICAL], view=[-10..50,-50..300]): print(p1);

<Text-field style="_cstyle277" layout="Bullet Item">Designing a Model</Text-field>

We are asked to compare a quadratic and a cubic model, so now we use the fit[leastsquare]( ) command to find the best-fit function for each model.

temp:= [[seq(data[i,1],i=1..nops(data))],[seq(data[i,2],i=1..nops(data))]]: fit[leastsquare[[x,y], y=a*x^2.,{a}]](temp): y2:=unapply(rhs(%),x);

fit[leastsquare[[x,y], y=a*x^3.,{a}]](temp): y3:=unapply(rhs(%),x);

p2:=plot(y2(x),x=0..45, labels=[xlabel, ylabel], color=blue, title="Quadratic Model", labeldirections=[HORIZONTAL, VERTICAL], view=[-10..50,-50..300]): print(p2);

p3:=plot(y3(x),x=0..45, labels=[xlabel, ylabel],title="Cubic Model", color=green, labeldirections=[HORIZONTAL, VERTICAL], view=[-10..50,-50..300]): print(p3);

print(display({p1,p2,p3}, labels=[xlabel, ylabel], title="Combined Model"));

The graphs seem to show that the cubic model better fits the data set .

<Text-field style="_cstyle278" layout="Bullet Item">Assessing the Errors</Text-field>

Now we analyze the fit for each of the models by calculating the percent errors, first for the quadratic function and then for the cubic.

<Text-field style="_cstyle293" layout="Heading 1">Analysis of the Errors for the Quadratic-Fit Function</Text-field>

Let's use the model to calculate the volume of lumber for each girth measurement in the original data set, and then calculate the residual errors and plot them.

quadpredictedvalues:=[seq([data[i,1], y2(data[i,1])], i=1..nops(data))];

quadresiduals:=[seq([data[i,1], data[i,2]-quadpredictedvalues[i,2]],i=1..nops(data))];

pointplot(quadresiduals, labels=[xlabel,"residual error"], labeldirections=[HORIZONTAL, VERTICAL]);

Note that the errors aren't completely random. There are more negative errors than there are positive ones, and the error seems to increase as the girth increases. These observations suggest that the quadratic function is not appropriate for modeling the relation between the volume of usable lumber in a tree and its girth.

It is usually more helpful to look at the error in relation to the size of the quantity being estimated. We calculate the relative percent errors as follows.

quadpercenterrors:=[seq([data[i,1],100*(data[i,2]-quadpredictedvalues[i,2])/data[i,2]],i=1..nops(data))];

Now we plot the percent errors and save the graph for a later comparison with the percent errors for the cubic function.

p3:=pointplot(quadpercenterrors, labels=[xlabel,"percent error"], labeldirections=[HORIZONTAL, VERTICAL],view=[-5..45,-150..20]): print(p3);

The percent error of largest magnitude is about -140%, and there is a trend in the errors, indicating that they are not random.

<Text-field style="_cstyle292" layout="Heading 1">Analysis of the Errors for the Cubic-Fit Function</Text-field>

Let's look at the errors for the cubic model. We do the same as we did for the quadratic function.

cubicpredictedvalues:=[seq([data[i,1], y3(data[i,1])], i=1..nops(data))];

cubicresiduals:=[seq([data[i,1], data[i,2]-cubicpredictedvalues[i,2]],i=1..nops(data))];

pointplot(cubicresiduals, labels=[xlabel,"residual error"], labeldirections=[HORIZONTAL, VERTICAL]);

Note that the residuals exhibit a more random pattern than they did in the quadratic model. There are as many negative errors as positive ones, and there is no apparent trend in the errors.

cubicpercenterrors:=[seq([data[i,1],100*(data[i,2]-cubicpredictedvalues[i,2])/data[i,2]],i=1..nops(data))];

p4:=pointplot(cubicpercenterrors, labels=[xlabel,"percent error"], color=blue, symbol=CROSS, labeldirections=[HORIZONTAL, VERTICAL]): print(p4);

In this case, the percent error of largest magnitude is about -20%, which is much better than for the quadratic regression function.

<Text-field style="_cstyle279" layout="Normal">Compare the Errors</Text-field>

To confirm our earlier assessment that the cubic regression is evidently better, we now plot the percent errors for the two models together.

print(display({p3,p4}, view=[-5..50,-150..40]));

<Text-field style="_cstyle294" layout="Bullet Item">Explaining the Model</Text-field>

The unit of board feet is a measure of volume. If a tree is modeled as a right circular cone, its volume is approximated by LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJWRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiPUYnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZRLUkmbWZyYWNHRiQ2KC1GIzYoLUklbXN1cEdGJDYlLUYsNiVRJSZwaTtGJy9GNUZCRj4tSSNtbkdGJDYkUSIyRidGPi8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRictRjs2LVExJkludmlzaWJsZVRpbWVzO0YnRj5GQEZDRkVGR0ZJRktGTS9GUFEmMC4wZW1GJy9GU0Zlby1GLDYlUSJyRidGNEY3RmFvLUYsNiVRImhGJ0Y0RjdGPi1GIzYkLUZbbzYkUSIzRidGPkY+LyUubGluZXRoaWNrbmVzc0dRIjFGJy8lK2Rlbm9tYWxpZ25HUSdjZW50ZXJGJy8lKW51bWFsaWduR0ZncC8lKWJldmVsbGVkR0ZCRj5GK0Y+ where V is the volume, h is the height of the tree, and r is its radius. The girth, g, is the circumference of the tree near the base so that LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Jy1GLDYlUSJnRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiPUYnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZRLUYjNigtSSNtbkdGJDYkUSIyRidGPi1GOzYtUTEmSW52aXNpYmxlVGltZXM7RidGPkZARkNGRUZHRklGS0ZNL0ZQUSYwLjBlbUYnL0ZTRmhuLUYsNiVRJSZwaTtGJy9GNUZCRj5GWi1GLDYlUSJyRidGNEY3Rj5GK0Y+RitGPg== or LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJyRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiPUYnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZRLUkmbWZyYWNHRiQ2KC1GIzYkLUYsNiVRImdGJ0Y0RjdGPi1GIzYmRistRiM2Ji1JI21uR0YkNiRRIjJGJ0Y+LUY7Ni1RMSZJbnZpc2libGVUaW1lcztGJ0Y+RkBGQ0ZFRkdGSUZLRk0vRlBRJjAuMGVtRicvRlNGYm8tRiw2JVElJnBpO0YnL0Y1RkJGPkY+RitGPi8lLmxpbmV0aGlja25lc3NHUSIxRicvJStkZW5vbWFsaWduR1EnY2VudGVyRicvJSludW1hbGlnbkdGXXAvJSliZXZlbGxlZEdGQkY+RitGPg== . If we assume that as a tree grows the proportion LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1JJm1mcmFjR0YkNigtRiM2JC1GLDYlUSJoRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnL0Y9USdub3JtYWxGJy1GIzYkLUYsNiVRInJGJ0Y5RjxGPy8lLmxpbmV0aGlja25lc3NHUSIxRicvJStkZW5vbWFsaWduR1EnY2VudGVyRicvJSludW1hbGlnbkdGSy8lKWJldmVsbGVkR1EmZmFsc2VGJy1JI21vR0YkNi1RIj1GJ0Y/LyUmZmVuY2VHRlAvJSpzZXBhcmF0b3JHRlAvJSlzdHJldGNoeUdGUC8lKnN5bW1ldHJpY0dGUC8lKGxhcmdlb3BHRlAvJS5tb3ZhYmxlbGltaXRzR0ZQLyUnYWNjZW50R0ZQLyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGX28tRiw2JVEia0YnRjlGPEY/RitGPw== is a constant, then we have that LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJWRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiPUYnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZRLUkmbWZyYWNHRiQ2KC1GIzYkLUYsNiVRJSZwaTtGJy9GNUZCRj5GPi1GIzYkLUkjbW5HRiQ2JFEiM0YnRj5GPi8lLmxpbmV0aGlja25lc3NHUSIxRicvJStkZW5vbWFsaWduR1EnY2VudGVyRicvJSludW1hbGlnbkdGYm8vJSliZXZlbGxlZEdGQkY+RitGPg== LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JC1JJW1zdXBHRiQ2JS1JKG1mZW5jZWRHRiQ2JC1JJm1mcmFjR0YkNigtRiM2JC1GLDYlUSJnRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnL0ZDUSdub3JtYWxGJy1GIzYmRistRiM2Ji1JI21uR0YkNiRRIjJGJ0ZFLUkjbW9HRiQ2LVExJkludmlzaWJsZVRpbWVzO0YnRkUvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRlUvJSlzdHJldGNoeUdGVS8lKnN5bW1ldHJpY0dGVS8lKGxhcmdlb3BHRlUvJS5tb3ZhYmxlbGltaXRzR0ZVLyUnYWNjZW50R0ZVLyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGXm8tRiw2JVElJnBpO0YnL0ZARlVGRUZFRitGRS8lLmxpbmV0aGlja25lc3NHUSIxRicvJStkZW5vbWFsaWduR1EnY2VudGVyRicvJSludW1hbGlnbkdGam8vJSliZXZlbGxlZEdGVUZFRksvJTFzdXBlcnNjcmlwdHNoaWZ0R1EiMEYnRkVGK0ZF LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkmbWZyYWNHRiQ2KC1GIzYkLUkjbWlHRiQ2JVEja2dGJy8lJ2l0YWxpY0dRJXRydWVGJy8lLG1hdGh2YXJpYW50R1EnaXRhbGljRicvRjhRJ25vcm1hbEYnLUYjNiYtRjE2I1EhRictRiM2Ji1JI21uR0YkNiRRIjJGJ0Y6LUkjbW9HRiQ2LVExJkludmlzaWJsZVRpbWVzO0YnRjovJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRk0vJSlzdHJldGNoeUdGTS8lKnN5bW1ldHJpY0dGTS8lKGxhcmdlb3BHRk0vJS5tb3ZhYmxlbGltaXRzR0ZNLyUnYWNjZW50R0ZNLyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGZm4tRjE2JVElJnBpO0YnL0Y1Rk1GOkY6Rj5GOi8lLmxpbmV0aGlja25lc3NHUSIxRicvJStkZW5vbWFsaWduR1EnY2VudGVyRicvJSludW1hbGlnbkdGYm8vJSliZXZlbGxlZEdGTUY6 = LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkmbWZyYWNHRiQ2KC1GIzYkLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GOFEnbm9ybWFsRictRiM2Ji1GMTYjUSFGJy1GIzYnLUkjbW5HRiQ2JFEjMjRGJ0Y6LUkjbW9HRiQ2LVExJkludmlzaWJsZVRpbWVzO0YnRjovJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRk0vJSlzdHJldGNoeUdGTS8lKnN5bW1ldHJpY0dGTS8lKGxhcmdlb3BHRk0vJS5tb3ZhYmxlbGltaXRzR0ZNLyUnYWNjZW50R0ZNLyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGZm4tRiM2JC1JJW1zdXBHRiQ2JS1GMTYlUSUmcGk7RicvRjVGTUY6LUZENiRRIjJGJ0Y6LyUxc3VwZXJzY3JpcHRzaGlmdEdRIjBGJ0Y6Rj5GOkY+RjovJS5saW5ldGhpY2tuZXNzR1EiMUYnLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRl1wLyUpYmV2ZWxsZWRHRk1GOg== LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JC1JJW1zdXBHRiQ2JS1GLDYlUSJnRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW5HRiQ2JFEiM0YnL0Y7USdub3JtYWxGJy8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGQUYrRkE= , where the last quantity LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkmbWZyYWNHRiQ2KC1GIzYkLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GOFEnbm9ybWFsRictRiM2Ji1GMTYjUSFGJy1GIzYnLUkjbW5HRiQ2JFEjMjRGJ0Y6LUkjbW9HRiQ2LVExJkludmlzaWJsZVRpbWVzO0YnRjovJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRk0vJSlzdHJldGNoeUdGTS8lKnN5bW1ldHJpY0dGTS8lKGxhcmdlb3BHRk0vJS5tb3ZhYmxlbGltaXRzR0ZNLyUnYWNjZW50R0ZNLyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGZm4tRiM2JC1JJW1zdXBHRiQ2JS1GMTYlUSUmcGk7RicvRjVGTUY6LUZENiRRIjJGJ0Y6LyUxc3VwZXJzY3JpcHRzaGlmdEdRIjBGJ0Y6Rj5GOkY+RjovJS5saW5ldGhpY2tuZXNzR1EiMUYnLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRl1wLyUpYmV2ZWxsZWRHRk1GOg== , is a constant. This shows that cubic relation between the volume of lumber produced and the girth of the tree near its base is a rational model.

<Text-field style="Heading 1" layout="Heading 1">Part V: Black Bass</Text-field>

<Text-field style="_cstyle280" layout="Bullet Item">Observing the Data</Text-field>

In the table that follows, L represents the lengths of New York black bass measured in inches (in.), and w represents the weight of the fish. Formulate and test a model that assumes the weight of the fish is proportional to the cube of its length. Can you provide a rational explanation of the model?

data := [[12.5, 17], [12.63, 16], [14.13, 17], [14.5, 23], [14.5, 26], [14.5, 27], [17.25, 41], [17.75, 49]];

matrix([[`L (in)`, `W (oz)`],op(data)]);

Next, we plot the data to see if there is a recognizable pattern.

xlabel:=`L(in)`: ylabel:=`w(oz)`: p1:=pointplot(data, labels=[xlabel, ylabel], view=[0..20,-5..50], labeldirections=[HORIZONTAL, VERTICAL]): print(p1);

<Text-field style="_cstyle281" layout="Bullet Item">Designing a Model</Text-field>

We are asked to construct a cubic model, so now we find the best-fit function using the fit[leastsquare]( ) command.

temp := [[seq(data[i,1],i=1..nops(data))],[seq(data[i,2],i=1..nops(data))]]: eqn:=fit[leastsquare[[x,y],y=a*x^3]]([temp[1],temp[2]]); y:=unapply(rhs(%),x):

Next, we plot the model function and superimpose its graph on the graph of the data.

p2:=plot(y(x), x=0..20, labels=[xlabel,ylabel], color=blue, labeldirections=[HORIZONTAL, VERTICAL], view=[0..20,-5..50]): print(p2);

print(display({p1,p2}));

The graphs show that the cubic model fits the data well.

<Text-field style="_cstyle282" layout="Bullet Item">Assessing the Errors</Text-field>

Now we analyze the fit by calculating the percent errors. First, we use the model to calculate the weight of the fish for each length measurement in the original data set, and then we calculate the residuals and plot them.

predictvalues:=[seq([data[i,1], y(data[i,1])], i=1..nops(data))];

res:=data-predictvalues: residuals:=[seq([data[i,1], res[i,2]],i=1..nops(data))];

pointplot(residuals, labels=[xlabel, `residual error`], labeldirections=[HORIZONTAL, VERTICAL]);

Note that there are nearly equal numbers of positive and negative errors, and there is no apparent trend in the errors. They appear to be random. As before, it is more meaningful to look at the error in relation to the size of the quantity being estimated, so we calculate the relative percent errors and plot them.

percenterrors:=[seq([data[i,1],100*(data[i,2]-predictvalues[i,2])/data[i,2]],i=1..nops(data))];

p3:=pointplot(percenterrors, labels=[xlabel,`% errors`], color=red, symbol=CROSS, labeldirections=[HORIZONTAL, VERTICAL]): print(p3);

<Text-field style="_cstyle283" layout="Bullet Item">Improving the Model</Text-field>

The largest percent error about -36%, but this data point seems to be an outlier. Let's see what happens if we remove this point from the data set.

data := [[12.5, 17], [12.63, 16], [14.13, 17], [14.5, 23], [14.5, 26], [14.5,27], [17.25, 41], [17.75, 49]]: data:=subsop(3=NULL,data); p1:=pointplot(data, view=[0..20,0..50], symbol=CROSS, color=blue,labels=[xlabel,ylabel], labeldirections=[HORIZONTAL, VERTICAL]): print(p1);

temp := [[seq(data[i,1],i=1..nops(data))],[seq(data[i,2],i=1..nops(data))]]: eqn:=fit[leastsquare[[x,y],y=a*x^3]]([temp[1],temp[2]]); y:=unapply(rhs(%),x): p2:=plot(y(x), x=0..20, labels=[xlabel,ylabel], color=red, labeldirections=[HORIZONTAL, VERTICAL]): print(p2);

print(display({p1,p2}));

predictvalues:=[seq([data[i,1], y(data[i,1])], i=1..nops(data))]: res:=data-predictvalues: residuals:=[seq([data[i,1], res[i,2]],i=1..nops(data))]: percenterrors:=[seq([data[i,1],100*(data[i,2]-predictvalues[i,2])/data[i,2]],i=1..nops(data))]; p3:=pointplot(percenterrors, labels=[xlabel,`% errors`], color=black, symbol=CROSS, labeldirections=[HORIZONTAL, VERTICAL]): print(p3);

After deleting the one data point, we find that the new model gives a largest percent error of about -11%, which is considerably better.

<Text-field style="_cstyle284" layout="Bullet Item">Explaining the Model</Text-field>

The weight of a fish is proportional to its volume. If we think of a fish to be shaped roughly like an ellipsoid, then its volume would be approximated by LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJWRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiPUYnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZRLUkmbWZyYWNHRiQ2KC1GIzYsLUkjbW5HRiQ2JFEiNEYnRj4tRjs2LVExJkludmlzaWJsZVRpbWVzO0YnRj5GQEZDRkVGR0ZJRktGTS9GUFEmMC4wZW1GJy9GU0Zbby1GLDYlUSUmcGk7RicvRjVGQkY+RmduLUYsNiVRIkxGJ0Y0RjdGZ24tRiw2JVEid0YnRjRGN0Znbi1GLDYlUSJkRidGNEY3Rj4tRiM2JC1GWjYkUSIzRidGPkY+LyUubGluZXRoaWNrbmVzc0dRIjFGJy8lK2Rlbm9tYWxpZ25HUSdjZW50ZXJGJy8lKW51bWFsaWduR0ZkcC8lKWJldmVsbGVkR0ZCRj5GK0Y+ where L, w, and d are its length, width, and depth, respectively. If, as the fish grows, we assume that its proportions remain the same then LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Jy1GLDYlUSJ3RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiPUYnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZRLUYjNiYtSSVtc3ViR0YkNiUtRiw2JVEia0YnRjRGNy1GIzYkRjFGPi8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnLUY7Ni1RMSZJbnZpc2libGVUaW1lcztGJ0Y+RkBGQ0ZFRkdGSUZLRk0vRlBRJjAuMGVtRicvRlNGX28tRiw2JVEiTEYnRjRGN0Y+RitGPkYrRj4= and LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Jy1GLDYlUSJkRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiPUYnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZRLUYjNiYtSSVtc3ViR0YkNiUtRiw2JVEia0YnRjRGNy1GIzYkRjFGPi8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnLUY7Ni1RMSZJbnZpc2libGVUaW1lcztGJ0Y+RkBGQ0ZFRkdGSUZLRk0vRlBRJjAuMGVtRicvRlNGX28tRiw2JVEiTEYnRjRGN0Y+RitGPkYrRj4=, where LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYkLUYvNiVRIkxGJ0YyRjUvRjZRJ25vcm1hbEYnLyUvc3Vic2NyaXB0c2hpZnRHUSIwRidGPQ== and LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYkLUYvNiVRImRGJ0YyRjUvRjZRJ25vcm1hbEYnLyUvc3Vic2NyaXB0c2hpZnRHUSIwRidGPQ== are constants. Therefore, V = LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1JJm1mcmFjR0YkNigtRiM2LC1JI21uR0YkNiRRIjRGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRictSSNtb0dGJDYtUTEmSW52aXNpYmxlVGltZXM7RidGOi8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQy8lKXN0cmV0Y2h5R0ZDLyUqc3ltbWV0cmljR0ZDLyUobGFyZ2VvcEdGQy8lLm1vdmFibGVsaW1pdHNHRkMvJSdhY2NlbnRHRkMvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZSLUYsNiVRJSZwaTtGJy8lJ2l0YWxpY0dGQ0Y6Rj0tRiw2JVEiTEYnL0ZZUSV0cnVlRicvRjtRJ2l0YWxpY0YnRj0tRiw2JVEmKGt3TClGJ0ZnbkZpbkY9LUYsNiVRJihrZEwpRidGZ25GaW5GOi1GIzYkLUY3NiRRIjNGJ0Y6RjovJS5saW5ldGhpY2tuZXNzR1EiMUYnLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRltwLyUpYmV2ZWxsZWRHRkMtRj42LVEiPUYnRjpGQUZERkZGSEZKRkxGTi9GUVEsMC4yNzc3Nzc4ZW1GJy9GVEZkcC1GMjYoLUYjNipGNkY9RlVGPS1GLDYlUSNrd0YnRmduRmluRj0tRiw2JVEja2RGJ0ZnbkZpbkY6RmFvRmZvRmlvRlxwRl5wRjpGK0Y6 LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JC1JJW1zdXBHRiQ2JS1GLDYlUSJMRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW5HRiQ2JFEiM0YnL0Y7USdub3JtYWxGJy8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGQUYrRkE= where the last quantity in parentheses is a constant. Therefore, the weight is proportional to the volume, which is proportional to the length cubed, and our model has a rational basis.

<Text-field style="Heading 1" layout="Heading 1">Part VI: Heart Rates of Mammals</Text-field>

<Text-field style="_cstyle261" layout="Bullet Item">Observing the Data</Text-field>

The following data relate the weight in grams (g) of some mammals to their heart rate in beats per minute. Plot the data. Is there a trend? If so, find a function that captures the trend of the data. Hint: try the form LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJ5RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiPUYnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZRLUklbXN1cEdGJDYlLUYsNiVRInhGJ0Y0RjctRiw2JVElLTEvbkYnRjRGNy8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGPkYrRj4= for n an integer.

data := [[4, 660], [25, 670], [200, 420], [300, 300], [2000, 205], [5000,120], [30000, 85], [50000, 70], [70000, 72], [450000, 38], [500000,40], [3000000, 48]];

matrix([[` x (g) `,`y (bpm) `],op(data)]);

Next, we plot the data to see if there is a recognizable pattern.

xlabel := "x(g)": ylabel := "y(bpm)": p1:=pointplot(data, symbol=CROSS, color=red, labels=[xlabel, ylabel], view=[-5000..70000,-50..700], labeldirections=[HORIZONTAL, VERTICAL]): print(p1);

Note that we did not plot the data points for the horse, the ox, and the elephant because that makes the scale too large to see if there is a pattern for the smaller mammals; however, these data points are included in the calculations that follow.

<Text-field style="_cstyle337" layout="Bullet Item">Designing a Model</Text-field>

The graph of the data suggests that there is a relationship. As suggested in the hint above, we now build a group of functions of the form LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJ5RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiPUYnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZRLUklbXN1cEdGJDYlLUYsNiVRInhGJ0Y0RjctRiw2JVElLTEvbkYnRjRGNy8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGPkYrRj4= where n is an integer. For integer values of n (1, 2, 3, 4, 5), we use the fit[leastsquare]( ) function to find the function of the form LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1cEdGJDYlLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GLzYlUSUtMS9uRidGMkY1LyUxc3VwZXJzY3JpcHRzaGlmdEdRIjBGJy9GNlEnbm9ybWFsRic= that best fits the data.

temp := [[seq(data[i,1],i=1..nops(data))],[seq(data[i,2],i=1..nops(data))]]: for n from 1 to 5 do eqn:=fit[leastsquare[[x,y],y=a*x^(-1.0/n),{a}]]([temp[1],temp[2]]): y[n]:=unapply(rhs(%),x); od;

p2:=plot([y[1](x), y[2](x),y[3](x),y[4](x),y[5](x)], x=0..70000, y=0..700, labels=[xlabel, ylabel],tickmarks=[4,7], color=[red,green,blue,black,brown], view=[-5000..70000,-50..700], labeldirections=[HORIZONTAL, VERTICAL]): print(p2);

Now we plot the models together with the data points.

print(display({p1,p2}));

<Text-field style="_cstyle342" layout="Bullet Item">Assessing the Errors</Text-field>

It appears that LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYkLUkjbW5HRiQ2JFEiNEYnL0Y2USdub3JtYWxGJ0Y+LyUvc3Vic2NyaXB0c2hpZnRHUSIwRidGPg== provides the best fit. Now we will analyze the errors. First, we calculate the heart rate values that our selected model predicts for each mammal, and then we calculate the residual errors and plot them.

predictedvalues:=[seq([data[i,1], y[4](data[i,1])], i=1..nops(data))];

residuals:=[seq([i, data[i,2]-predictedvalues[i,2]],i=1..nops(data))];

pointplot(residuals, color=red, labels=[xlabel,`residual error`], view=[-2..12,-175..175], labeldirections=[HORIZONTAL, VERTICAL]);

As before, it is helpful to look at the error in relation to the size of the quantity being estimated. We calculate the relative percent errors.

percenterrors:=[seq([i,100*(data[i,2]-predictedvalues[i,2])/data[i,2]],i=1..nops(data))];

pointplot(percenterrors, color=red, labels=[xlabel,`percent error`], labeldirections=[HORIZONTAL, VERTICAL]);

The errors appear to be somewhat random, but there is possibly a trend in the errors, with two outliers. Can you see the trend and the outliers? The largest relative percentage errors are for the largest and the smallest animals (i.e., the bat and the elephant) with magnitudes of 23% and 42%, respectively. The model appears to capture a trend in the data, which could be useful in understanding the relationship between mammal size and heart rate; however, it probably would not be useful as a predictive tool since the magnitudes of the individual errors are so large.

<Text-field style="Heading 1" layout="Heading 1">You Try It: A Rational Explanation for Heart Rates</Text-field>

Try to provide a rational explanation for the heart rate model that we found in Part VI.

<Text-field style="Heading 1" layout="Heading 1">You Try It: Modeling Takes Practice</Text-field>

Mathematical modeling and a computer algebra system like Maple help us describe real-world phenomena, understand the mechanisms behind the behavior, and make predictions. The procedure for modeling generally includes the following steps:

1. Plot and observe the data, looking for relationships and patterns. 2. Formulate a function to model the data. 3. Assess your model by analyzing the residual errors and/or relative errors. 4. Improve your model, if possible. 5. Use your model to gain a better understanding of the phenomenon you are modeling.

6. Use your model to make predictions.

With some practice, you too can become an effective mathematical modeler. There are plenty of modeling exercises that you can do yourself. Select a problem and obtain some data (possibly by designing your own experiment) and use the modeling procedures outlined in the preceding Parts of this module as a template for mathematical modeling. To help you get started, we include the data sets for four modeling examples.

Drug Levels

drugdata := [[0, 853], [1, 587], [2, 390], [3, 274], [4, 189], [5, 130], [6,97], [7, 67], [8, 50], [9, 40], [10, 31]];

Cell Count

celldata:= [[0, 597], [2, 893], [4, 1339], [6, 1995], [8, 2976], [10, 4433], [12, 6612], [14, 9865], [16, 14719], [18, 21956], [20, 32763]];

Vacuum Pump

vacuumdata := [[0, 100000], [1, 36788], [2, 13537], [3, 4986], [4, 1837], [5, 671]];

Doctoral Degrees

doctoraldata := [[6, 520], [10, 460], [14, 680], [18, 630], [20, 730], [21,810], [22, 830]];

LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic=