/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 11.08.0 ] */

/* [wxMaxima: title   start ]
Modeling Change: Springs, Driving 
Safety, Radioactivity, Trees, Fish, and
Mammals
   [wxMaxima: title   end   ] */

/* [wxMaxima: comment start ]
(Chapter 1 - Technology Application Projects)
   [wxMaxima: comment end   ] */

/* [wxMaxima: section start ]
Introduction
   [wxMaxima: section end   ] */

/* [wxMaxima: comment start ]
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 Maxima is a powerful tool that will aid you 
in your mathematical modeling.

In this worksheet, we use the two Maxima packages draw and 
stats. We initialize them at the beginning Part I.
   [wxMaxima: comment end   ] */

/* [wxMaxima: subsect start ]
Technology Guidelines
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
NOTE:  If you have just finished a document, restart Maxima 
before executing a new document.  This can be done by choosing
"Restart" from the Maxima menu.

TO OPEN OR CLOSE CELLS
Click on the arrow at the top of the cell bracket. 

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

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

Alternatively, you can execute the entire worksheet by 
selecting the "Evaluate All Cells" command from the "Cell" drop
down menu or simply press Ctrl-r. 

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 Maxima and start it up 
again. 
   [wxMaxima: comment end   ] */

/* [wxMaxima: section start ]
Part I:  Spring Elongation 
   [wxMaxima: section end   ] */

/* [wxMaxima: subsect start ]
Observing the Data
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
load(draw)$
load(stats)$
load(lsquares)$
ratprint: false$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
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]];
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
datamatrix: matrix(["mass","stretch (in)"],[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]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Next, we plot the data to see if there is a recognizable pattern
and name the plot p1 for later use.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
p1: wxdraw2d(
    nticks=200,
    color=blue,
    point_type=6,
    point_size=1,
    points(data),
    xlabel="mass",
    ylabel="stretch (in)");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
Designing a Model
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
The data strongly suggest a linear relationship. Now we use 
the simple_linear_regression( ) function to find the line of 
best-fit or the linear regression line.

The output of the function simple_linear_regression is an 
inference_result Maxima object with the following results: 

1.  model: the equation for the line of best fit.   
2.  means: bivariate mean. 
3.  variances: variances of both variables. 
4.  correlation: correlation coefficient. 
5.  adc: adjusted determination coefficient. 
6.  a_estimation: estimation of parameter a. 
7.  a_conf_int: confidence interval of parameter a. 
8.  b_estimation: estimation of parameter b. 
9.  b_conf_int: confidence interval of parameter b.
10. hypotheses: null and alternative hypotheses about parameter b. 
11. statistic: value of the sample statistic used for testing the 
    null hypothesis. 
12. distribution: distribution of the sample statistic, together 
    with its parameter. 
13. p_value: p-value of the test about b. 
14. v_estimation: unbiased variance estimation, or residual 
    variance. 
15. v_conf_int: variance confidence interval. 
16. cond_mean_conf_int: confidence interval for the conditioned 
    mean.
17. new_pred_conf_int: confidence interval for a new prediction.
18. residuals: list of pairs (prediction, residual), ordered 
    with respect to predictions. This is useful for goodness of 
    fit analysis.

Only items 1, 4, 14, 9, 10, 11, 12, and 13 above, in this order,
are shown by default. The rest remain hidden until the user 
makes use of functions items_inference and take_inference. 
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
data2: simple_linear_regression(data);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
eqn: take_inference(model,data2);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
The constant of proportionality is 0.8742. 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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
p2: wxdraw2d(
    nticks=200,
    color=red,
    explicit(eqn,x,0,10),
    xlabel="mass",
    ylabel="stretch (in)");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    color=blue,
    point_type=6,
    point_size=1,
    points(data),
    color=red,
    explicit(eqn,x,0,10),
    xlabel="mass",
    ylabel="stretch (in)");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
Assesssing the Errors
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
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.

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

Note that a list of predicted value / residual pairs is
automatically created when the simple_linear_regression( ) 
function is executed.  So we will retrieve this list, and then
save just the residual portion as a new list.  Next we create
the residual plot.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
kill(residuals);
PredictValues_Residuals: take_inference(residuals,data2);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
residuals: makelist([i,PredictValues_Residuals[i+1][2]],i,0,10);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    yrange=[-.04,.04],
    xrange=[-1,11],
    xlabel="mass",
    ylabel="residual error", 
    color=blue,
    point_type=6,
    point_size=1,
    points(residuals));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
percenterrors: makelist([i,100*residuals[i+1][2]/data[i+1][2]],i,1,10);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
percenterrplot: points(percenterrors);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    yrange=[-2,1],
    xrange=[-2,12],
    xlabel="mass",
    ylabel="% error", 
    color=blue, percenterrplot);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: section start ]
Part II: Safe Following Distance
   [wxMaxima: section end   ] */

/* [wxMaxima: subsect start ]
Observing the Data
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
kill(all)$
load(draw)$
load(stats)$
load(lsquares)$
ratprint: false$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
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]];
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
datalabelsmatrix: matrix(["speed","stopping dist"],[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]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
datamatrix: matrix([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]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
We plot the data to see if there is a recognizable pattern.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
p1: wxdraw2d(
    xlabel="speed (mph)",
    ylabel="stopping distance (ft)",
    color=blue,
    point_type=6,
    point_size=1,
    points(data));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
Designing a Model
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
The data suggest a possible quadratic relationship. Now we load 
the lsquares package in order to run a nonlinear regression.  
We will use the lsquares_estimates function to find the curve of 
best-fit.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
load(lsquares)$
numer:true$
ratprint: false$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
coeff: lsquares_estimates(datamatrix, [x,y], y=a+b*x+c*x^2, [a,b,c]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
a: rhs(coeff[1][1]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
b: rhs(coeff[1][2]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
c: rhs(coeff[1][3]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
eqn: a+b*x+c*x^2;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Next, we plot the regression curve and then show the plot of the
data and the model together on the same graph.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
p2: wxdraw2d(
    nticks=10,
    color=red,
    explicit(eqn,x,0,85),
    xlabel="speed (mph)",
    ylabel="stopping distance (ft)");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=10,
    color=red,
    explicit(eqn,x,0,85),
    color=blue,
    point_type=6,
    point_size=1,
    points(data),
    xlabel="speed (mph)",
    ylabel="stopping distance (ft)");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
The best-fit quadratic function appears to fit the data very 
well.
   [wxMaxima: comment end   ] */

/* [wxMaxima: subsect start ]
Assessing the Errors
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
predictvalues: makelist([i*5,subst(i*5,x,eqn)],i,4,15);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
predictvalueslabels: matrix(["speed","predicted values"],
    [20,54.47252747252747],[25,74.9000999000999],[30,98.55544455544455],
    [35,125.4385614385614],[40,155.5494505494505],[45,188.8881118881119],
    [50,225.4545454545454],[55,265.2487512487512],[60,308.2707292707293],
    [65,354.5204795204795],[70,403.998001998002],[75,456.7032967032967]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
  res: data-predictvalues;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
residuals: makelist([20+5*(i-1),res[i][2]],i,1,length(data));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    color=blue,
    point_type=6,
    point_size=1,
    points(residuals),
    xlabel="residual error");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
percenterrors: makelist([20+5*(i-1),100*residuals[i][2]/data[i][2]],i,2,11);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    color=blue,
    point_type=6,
    point_size=1,
    points(percenterrors),
    xlabel="percentage error");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
The maximum relatvie error of this model is about 0.6%.
   [wxMaxima: comment end   ] */

/* [wxMaxima: subsect start ]
Improving the Model
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
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.)
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
data2: makelist([data[i][1],(data[i][1]*88/60.0)*2,data[i][2]],i,1,12);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
matrix(["v (mph)","2-sec dist. (ft)","stopping dist (ft)"],[20,58.66666666666666,54],
    [25,73.33333333333333,75],[30,88.0,98],[35,102.6666666666667,126],
    [40,117.3333333333333,156],[45,132.0,190],[50,146.6666666666667,226],
    [55,161.3333333333333,265],[60,176.0,307],[65,190.6666666666667,354],
    [70,205.3333333333333,402],[75,220.0,459]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
data3: makelist([data[i][1], data[i][2]/(data[i][1]*88/60)],i,1,length(data));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
matrix(["speed (mph)","separation time (sec)"],[20,1.840909090909091],
    [25,2.045454545454545],[30,2.227272727272727],[35,2.454545454545455],
    [40,2.659090909090909],[45,2.878787878787879],[50,3.081818181818182],
    [55,3.285123966942149],[60,3.488636363636364],[65,3.713286713286713],
    [70,3.915584415584415],[75,4.172727272727273]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Now we plot the separation time versus speed.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    color=blue,
    point_type=6,
    point_size=1,
    points(data3),
    xlabel="speed (mph)",
    ylabel="separation time (sec)",
    xrange=[-5,80],
    yrange=[-1,5]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
linreg: simple_linear_regression(data3);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
eqn: take_inference(model,linreg);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    color=red,
    explicit(eqn,x,0,85),
    xrange=[-5,80],
    yrange=[-1,5],
    xlabel="speed (mph)",
    ylabel="separation time (sec)");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    color=blue,
    point_type=6,
    point_size=1,
    points(data3),
    color=red,
    explicit(eqn,x,0,85),
    xrange=[-5,80],
    yrange=[-1,5],
    xlabel="speed (mph)",
    ylabel="separation time (sec)");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
You Try It:  Make a Better Rules of Thumb

Using the Model
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: section start ]
Part III: Radioactivity
   [wxMaxima: section end   ] */

/* [wxMaxima: subsect start ]
Observing the Data
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
kill(all)$
load(draw)$
load(stats)$
load(lsquares)$
ratprint: false$
/* [wxMaxima: input   end   ] */

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

/* [wxMaxima: input   start ] */
datamatrix: matrix([0, 10023], [1, 8174], [2, 6693], [3, 5500], [4, 4489], 
    [5, 3683], [6,3061], [7, 2479], [8, 2045], [9, 1645], [10, 1326]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
matrix(["Time (min)", "Radioactivity (cpm)"],[0, 10023], [1, 8174], 
    [2, 6693], [3, 5500], [4, 4489], [5, 3683], [6,3061], [7, 2479], 
    [8, 2045], [9, 1645], [10, 1326]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
We begin by plotting the data points.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    color=blue,
    point_type=6,
    point_size=1,
    points(data),
    xlabel="x (min)",
    ylabel="y (cpm)",
    xrange=[-2,10],
    yrange=[-1000,11000]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
Designing a Model
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
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 lsquares_estimates function. What does the function 
look like to you?  A decaying exponential? We try to find a 
function of the form y=ae^(-kx), where we vary the value of k 
and the computer finds the a to find the best-fit function of 
this form.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
numer:true$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
coeff: lsquares_estimates(datamatrix, [x,y], y=a*%e^(-0.1*x), [a]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
a: rhs(coeff[1]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
f2: a*%e^(-0.1*x);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    color=red,
    explicit(f2,x,0,12),
    xlabel="x (min)",
    ylabel="y (cpm)",
    xrange=[-2,10],
    yrange=[-1000,11000]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    color=blue,
    point_type=6,
    point_size=1,
    points(data),
    color=red,
    explicit(f2,x,0,12),
    xlabel="x (min)",
    ylabel="y (cpm)",
    xrange=[-2,10],
    yrange=[-1000,11000]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
Assesssing the Errors
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
predictvalues: makelist([i,subst(i,x,f2)],i,0,10);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
res: data-predictvalues;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
residuals: makelist([i-1,res[i][2]],i,1,11);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
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) of all of the individual or local residuals, but 
this can be very misleading.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
residavg: sum(residuals[i][2],i,1,11)/length(residuals);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
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 
"least squares." Squaring the residuals removes the canceling 
effect that occurs when we add them together for a measure of 
the global error. The Maxima lsquares_estimates( ) 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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
residsquaresum: sum(residuals[i][2]^2,i,1,11);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Before we compare this value, we calculate the mean of the 
squared residuals.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
msresiduals: residsquaresum/length(residuals);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
rmsresiduals: sqrt(msresiduals);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
The value that we calculate in the preceding step is the root 
of the mean of the squares of the residuals and is often times 
called the root-mean-square or "rms" value of the residuals.
   [wxMaxima: comment end   ] */

/* [wxMaxima: subsect start ]
You Try It: Improving the Model and 
Making Predictions
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
kill(a);
k: 0.2;
coeff: lsquares_estimates(datamatrix, [x,y], y=a*%e^(-k*x), [a]);
a: rhs(coeff[1]);
f2: a*%e^(-k*x);
wxdraw2d(
    nticks=200,
    color=red,
    explicit(f2,x,0,10),
    color=blue,
    point_type=6,
    point_size=1,
    points(data),
    xlabel="x (min)",
    ylabel="y (cpm)",
    xrange=[-2,10],
    yrange=[-1000,11000]);
residuals: makelist([data[i][2],data[i][2]-subst(data[i][1],x,f2)],i,1,length(data));
residsquaresum: sum(residuals[i][2]^2,i,1,11);
mresiduals: residsquaresum/length(data);
rmsresiduals: sqrt(mresiduals);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
MAKING PREDICTIONS

Now use your improved model to determine when the radioactivity
will fall below 500 counts per minute.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
solve([f2=500], [x]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Your improved model should predict that the radioactivity will 
be below 500 cpm after about 15 minutes. Is that what you get?
   [wxMaxima: comment end   ] */

/* [wxMaxima: section start ]
Part IV: Ponderosa Pines
   [wxMaxima: section end   ] */

/* [wxMaxima: subsect start ]
Observing the Data
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
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?
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
kill(all)$
load(draw)$
load(stats)$
load(lsquares)$
ratprint: false$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
data: [[17, 19], [19, 25], [20, 32], [23, 57], [25, 71], [28, 113], [32,123], [38, 252], 
[39, 259], [41, 294]];
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
datamatrix: matrix([17, 19], [19, 25], [20, 32], [23, 57], [25, 71], [28, 113], [32,123], [38, 252], 
[39, 259], [41, 294]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
matrix(["girth (in)","lumber (bf)"],[17, 19], [19, 25], [20, 32], [23, 57], [25, 71], 
    [28, 113], [32,123], [38, 252], [39, 259], [41, 294]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    color=blue,
    point_type=6,
    point_size=1,
    points(data),
    xlabel="girth (in)",
    ylabel="lumber (bf)",
    xrange=[-10,50],
    yrange=[-50,300]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
Designing a Model
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
We are asked to compare a quadratic and a cubic model, so now 
we use the lsqaures_estimates( ) function to find the best-fit 
function for each model.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
coeff1: lsquares_estimates(datamatrix, [x,y], y=a*x^2, [a]);
a: rhs(coeff1[1]);
y2: a*x^2;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
coeff2: lsquares_estimates(datamatrix, [x,y], y=b*x^3, [b]);
b: rhs(coeff2[1]);
y3: b*x^3;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    title="Quadratic Model",
    color=blue,
    explicit(y2,x,0,45),
    xlabel="girth (in)",
    ylabel="lumber (bf)",
    xrange=[-10,50],
    yrange=[-50,300]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    title="Cubic Model",
    color=green,
    explicit(y3,x,0,45),
    xlabel="girth (in)",
    ylabel="lumber (bf)",
    xrange=[-10,50],
    yrange=[-50,300]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    title="Combined Model",
    color=blue,
    explicit(y2,x,0,45),
    color=green,
    explicit(y3,x,0,45),
    color=black,
    point_type=6,
    point_size=1,
    points(data),
    xlabel="girth (in)",
    ylabel="lumber (bf)",
    xrange=[-10,50],
    yrange=[-50,300]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
The graphs seem to show that the cubic model better fits the 
data set.
   [wxMaxima: comment end   ] */

/* [wxMaxima: subsect start ]
Assessing the Errors
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
ANALYSIS OF THE ERRORS FOR THE QUADRATIC-FIT FUNCTION

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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
quadpredictedvalues: makelist([data[i][1], subst(data[i][1],x,y2)],i,1,length(data));
quadresiduals: makelist([data[i][1],data[i][2]-quadpredictedvalues[i][2]],i,1,length(data));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    color=red,
    point_type=6,
    point_size=1,
    points(quadresiduals),
    xlabel="girth (in)",
    ylabel="residual error");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
quadpercenterrors: makelist([data[i][1],
    100*(data[i][2]-quadpredictedvalues[i][2])/data[i][2]],i,1,length(data));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Now we plot the percent errors.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    color=red,
    point_type=6,
    point_size=1,
    points(quadpercenterrors),
    xrange=[-5,45],
    yrange=[-150,20],
    xlabel="girth (in)",
    ylabel="percent error");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
The percent error of largest magnitude is about -140%, and there
is a trend in the errors, indicating that they are not random.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
ANALYSIS OF THE ERRORS FOR THE CUBIC-FIT FUNCTION

Let's look at the errors for the cubic model.  We do the same 
as we did for the quadratic function.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
cubicpredictedvalues: makelist([data[i][1], subst(data[i][1],x,y3)],i,1,length(data));
cubicresiduals: makelist([data[i][1],data[i][2]-cubicpredictedvalues[i][2]],i,1,length(data));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    color=blue,
    point_type=6,
    point_size=1,
    points(cubicresiduals),
    xlabel="girth (in)",
    ylabel="residual error");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
cubicpercenterrors: makelist([data[i][1],
    100*(data[i][2]-cubicpredictedvalues[i][2])/data[i][2]],i,1,length(data));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    color=blue,
    point_type=6,
    point_size=1,
    points(cubicpercenterrors),
    xlabel="girth (in)",
    ylabel="percent error");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
In this case, the percent error of largest magnitude is about 
-20%, which is much better than for the quadratic regression 
function.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
COMPARE THE ERRORS

To confirm our earlier assessment that the cubic regression is 
evidently better, we now plot the percent errors for the two 
models together.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    color=blue,
    point_type=6,
    point_size=1,
    points(cubicpercenterrors),
    color=red,
    points(quadpercenterrors),
    xlabel="girth (in)",
    ylabel="percent error");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
Explaining the Model
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
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
V = (pi*r^2*h)/3 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 g = 2*pi*r or r = g/(2*pi). If 
we assume that as a tree grows the proportion h/r = k is a 
constant, then we have that 
V = (pi/3)*(g/(2*pi))^2*(kg/(2*pi) = k/(24*pi^2)*g^3, where the 
last quantity k/(24*pi^2), 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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: section start ]
Part V: Black Bass
   [wxMaxima: section end   ] */

/* [wxMaxima: subsect start ]
Observing the Data
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
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?
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
kill(all)$
load(draw)$
load(stats)$
load(lsquares)$
ratprint: false$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
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]];
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
datamatrix: matrix([12.5, 17], [12.63, 16], [14.13, 17], [14.5, 23], [14.5, 26], [14.5, 27], 
    [17.25, 41], [17.75, 49]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
matrix(["L (in)", "W (oz)"],[12.5, 17], [12.63, 16], [14.13, 17], [14.5, 23], [14.5, 26], 
    [14.5, 27], [17.25, 41], [17.75, 49]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Next, we plot the data to see if there is a recognizable 
pattern.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    color=blue,
    point_type=6,
    point_size=1,
    points(data),
    xlabel="L (in)",
    ylabel="w (oz)",
    xrange=[0,20],
    yrange=[-5,50]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
Designing a Model
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
We are asked to construct a cubic model, so now we find the 
best-fit function using the lsquares_estimates( ) command.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
coeff: lsquares_estimates(datamatrix, [x,y], y=a*x^3, [a]);
a: rhs(coeff[1]);
y4: a*x^3;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Next, we plot the model function and superimpose its graph on 
the graph of the data.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    color=blue,
    explicit(y4,x,0,45),
    color=black,
    point_type=6,
    point_size=1,
    points(data),
    xlabel="L (in)",
    ylabel="w (oz)",
    xrange=[0,20],
    yrange=[-5,50]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
The graphs show that the cubic model fits the data well.
   [wxMaxima: comment end   ] */

/* [wxMaxima: subsect start ]
Assessing the Errors
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
cubicpredictedvalues: makelist([data[i][1], subst(data[i][1],x,y4)],i,1,length(data));
cubicresiduals: makelist([data[i][1],data[i][2]-cubicpredictedvalues[i][2]],i,1,length(data));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    color=blue,
    point_type=6,
    point_size=1,
    points(cubicresiduals),
    xlabel="L (in)",
    ylabel="residual error");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
cubicpercenterrors: makelist([data[i][1],
    100*(data[i][2]-cubicpredictedvalues[i][2])/data[i][2]],i,1,length(data));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    color=blue,
    point_type=6,
    point_size=1,
    points(cubicpercenterrors),
    xlabel="L (in)",
    ylabel="percent error");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
Improving the Model
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
data: [[12.5, 17], [12.63, 16], [14.5, 23], [14.5, 26], [14.5,27], 
    [17.25, 41], [17.75, 49]];
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
datamatrix: matrix([12.5, 17], [12.63, 16], [14.5, 23], [14.5, 26], [14.5,27], 
    [17.25, 41], [17.75, 49]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    color=blue,
    point_type=6,
    point_size=1,
    points(data),
    xrange=[0,20],
    yrange=[0,50],
    xlabel="L (in)",
    ylabel="w (oz)");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
kill(coeff,a);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
coeff: lsquares_estimates(datamatrix, [x,y], y=a*x^3, [a]);
a: rhs(coeff[1]);
y5: a*x^3;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    color=blue,
    explicit(y5,x,0,45),
    color=black,
    point_type=6,
    point_size=1,
    points(data),
    xlabel="L (in)",
    ylabel="w (oz)",
    xrange=[0,20],
    yrange=[-5,50]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
cubicpredictedvalues: makelist([data[i][1], subst(data[i][1],x,y5)],i,1,length(data));
cubicresiduals: makelist([data[i][1],data[i][2]-cubicpredictedvalues[i][2]],i,1,length(data));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
cubicpercenterrors: makelist([data[i][1],
    100*(data[i][2]-cubicpredictedvalues[i][2])/data[i][2]],i,1,length(data));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    color=blue,
    point_type=6,
    point_size=1,
    points(cubicpercenterrors),
    xlabel="L (in)",
    ylabel="percent error");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
After deleting the one data point, we find that the new model 
gives a largest percent error of about -11%, which is 
considerably better.
   [wxMaxima: comment end   ] */

/* [wxMaxima: subsect start ]
Explaining the Model
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
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 V = (4*pi*L*w*d)/3 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 w = k(w)*L and d = k(d)*L, where k(w) and k(d) are 
constants. Therefore, 
V = (4*pi*L*(k(w)*L)(k(d)*L)/3 = (4*pi*k(w)*k(d)/3)*L^3
where (4*pi*k(w)*k(d)/3) 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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: section start ]
Part VI: Heart Rates of Mammals
   [wxMaxima: section end   ] */

/* [wxMaxima: subsect start ]
Observing the Data
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
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 y = x^(-1/n) for n an integer.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
kill(all)$
load(draw)$
load(stats)$
load(lsquares)$
ratprint: false$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
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]];
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
datamatrix: matrix([4, 660], [25, 670], [200, 420], [300, 300], [2000, 205], 
    [5000,120], [30000, 85], [50000, 70], [70000, 72], [450000, 38], 
    [500000,40], [3000000, 48]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
matrix(["x (g)", "y (bpm)"],[4, 660], [25, 670], [200, 420], [300, 300], [2000, 205], 
    [5000,120], [30000, 85], [50000, 70], [70000, 72], [450000, 38], 
    [500000,40], [3000000, 48]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Next, we plot the data to see if there is a recognizable 
pattern.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    color=blue,
    point_type=6,
    point_size=1,
    points(data),
    xlabel="x (g)",
    ylabel="y (bpm)",
    xrange=[-5000,70000],
    yrange=[-50,700]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: subsect start ]
Designing a Model
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
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 y = x^(-1/n) where n is an integer. For integer 
values of n (1, 2, 3, 4, 5), we use the lsquares_estimates( ) 
function to find the function of the form x^(-1/n) that best 
fits the data.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
for n: 1 thru 5 do block(
    kill(a),
    coeff: lsquares_estimates(datamatrix, [x,y], y=a*x^(-1.0/n), [a]),
    a: rhs(coeff[1]),
    y[n]: a*x^(-1.0/n));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
y[1];
y[2];
y[3];
y[4];
y[5];
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    color=red,
    explicit(y[1],x,0,70000),
    color=green,
    explicit(y[2],x,0,70000),
    color=blue,
    explicit(y[3],x,0,70000),
    color=black,
    explicit(y[4],x,0,70000),
    color=brown,
    explicit(y[5],x,0,70000),
    xlabel="x (g)",
    ylabel="y (bpm)",
    xrange=[-5000,70000],
    yrange=[-50,700]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Now we plot the models together with the data points.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=200,
    color=red,
    explicit(y[1],x,0,70000),
    color=green,
    explicit(y[2],x,0,70000),
    color=blue,
    explicit(y[3],x,0,70000),
    color=black,
    explicit(y[4],x,0,70000),
    color=brown,
    explicit(y[5],x,0,70000),
    color=purple,
    point_type=6,
    point_size=1,
    points(data),
    xlabel="x (g)",
    ylabel="y (bpm)",
    xrange=[-5000,70000],
    yrange=[-50,700]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
Assessing the Errors
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
It appears that y[4] 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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
predictedvalues: makelist([data[i][1], subst(data[i][1],x,y[4])],i,1,length(data));
residuals: makelist([i,data[i][2]-predictedvalues[i][2]],i,1,length(data));
wxdraw2d(
    nticks=200,
    color=blue,
    point_type=6,
    point_size=1,
    points(residuals),
    xrange=[2,12],
    yrange=[-175,175],
    xlabel="x (g)",
    ylabel="residuals");
percenterrors: makelist([i,100*(data[i][2]-predictedvalues[i][2])/data[i][2]],i,1,length(data));
wxdraw2d(
    nticks=200,
    color=blue,
    point_type=6,
    point_size=1,
    points(percenterrors),
    xrange=[2,12],
    yrange=[-30,50],
    xlabel="x (g)",
    ylabel="percent error");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: subsect start ]
You Try It:  A Rational Explanation for Heart Rate
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
Try to provide a rational explanation for the heart rate model
that we created in Part VI.
   [wxMaxima: comment end   ] */

/* [wxMaxima: section start ]
You Try It:  Modeling Takes Practice
   [wxMaxima: section end   ] */

/* [wxMaxima: comment start ]
Mathematical modeling and a computer algebra system like Maxima
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.
   [wxMaxima: comment end   ] */

/* [wxMaxima: subsect start ]
Drug Levels
   [wxMaxima: subsect end   ] */

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

/* [wxMaxima: subsect start ]
Cell Count
   [wxMaxima: subsect end   ] */

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

/* [wxMaxima: subsect start ]
Vacuum Pump
   [wxMaxima: subsect end   ] */

/* [wxMaxima: input   start ] */
vacuumdata: [[0, 100000], [1, 36788], [2, 13537], [3, 4986], [4, 1837], 
    [5, 671]];
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
Doctoral Degrees
   [wxMaxima: subsect end   ] */

/* [wxMaxima: input   start ] */
doctoraldata: [[6, 520], [10, 460], [14, 680], [18, 630], [20, 730], 
    [21,810], [22, 830]];
/* [wxMaxima: input   end   ] */

/* Maxima can't load/batch files which end with a comment! */
"Created with wxMaxima"$