Start R in Calculus

DANIEL KAPLAN
S TA R T R IN CALCULUS
PROJECT MOSAIC
2
Copyright (c) 2013 by Daniel Kaplan.

Beta Edition, January 2013
All rights reserved. No part of the material protected
by this copyright notice may be reproduced or uti-
lized in any form, electronic or mechanical, including
photocopying, recording, or scanning, or stored on
any information storage or retrieval system, without
written permission from the author.
Cover Photo: . [Photo credit: Maya Hanna.]

3
To Netta, Liat, and Tamar.

Preface: Calculus and Calculation
Its no coincidence that calculus and calculate are

similar. When you perform a calculation, you take infor-
mation in one form say 3 19 and translate it into
another, more readily interpreted form 57. Similarly,
calculus translates information about parts and about
change into information about wholes and accumulation.
In the 17thINFINIMENT
DES PITITS.
and 18th centuries, /. Partle.
when calculus was in-
de la feconde ^.par la premiere xy ce qui donne x^
vented by Newton and others, todays familiar comput-
(
& partant
ers were la difference de xyz^ fera yz^dx -t- x
unimaginable. Sensibly, the calculations were
+-
presented
xydz^ in a means suitable to the technology of those
3. Lapaper
centuries:
difference de xyz.it eft uyzjlx -+- ux^dy
and quill. The calculus computations of
Ce qui fe prouve comme dans Ic
that era involved the
xygjiu.manipulation of symbols according
regardant le produit xyz^ comme
cas precedent en unc
tofeumathematicallyenderived rules, for instance x2 ! 2x.
eft ainfi des autres a rinfini d oul on
lequantite.il ,
forme cette regie.

RE G L E II. Figure 1: One of the many
rules for the manipulation
Pour la quantites multiplies}. of symbols presented in the
du produit de first calculus textbook in 1699,
La. difference plufieurs quantites multi-
by the Marquis de lHopital.
plieesles unes paries autres, eft egale alafommedespro- ("Regle" is French for "rule.")
duits de la difference de chacune de ces quantites par le
produit des autres.

Ainfi la difference de ax eft x o + adx, c eft a dire
adx. Celle de a -t- x x b

y&bdx ydx ady
xdy.
PROPOSITION III.
These rules are still useful in important ways, but there

Probleme.
are
6-now
IRENDRE/^ difference dperforming
other technologies for une frattion the computa-
quelconqus,
tions of calculus. In addition to symbolic calculus, there
La difference de eft yd ~/ Jy Car fuppofant j-
*
is numerical calculus,^ which relies on simple arithmetic.
=
=
.
<,
on aura xcalculus, 6c comme ces deux

Numerical y z^, tedious and impractically slow varia
quantites for
bles x&y z^ doivent toujours etre entr elles, fbic
a human, is easy work for a computer. Perhaps
egales surpris-
qu elles ou
ingly, numerical calculus isdiminuent,
augmentent il s enfuit
also easier for many people que leurto
difference, c eft d dire leurs
accroiffemensou diminutions
feront auili egales entr elles 5
partant &
* on aura dx
=yd^ -t-
*dy ^.d^=^~^LL
%
l^L=^ en mettanc
^ fa valeur ~ Ce il falloit &c. d ou Ton forme
pour qu
.
,
5
understand and carry out. As such, it provides a means

to teach the concepts of calculus.
Numerical calculus is also more general than symbolic
calculus. There are many sorts of problems for which
symbolic calculus fails to provide an answer even for the
most talented and well-trained mathematicians. By using
numerical calculus, a much wider range of problems can
be addressed, even by beginners. This is especially true
for modern problems in modeling and data analysis for
which the symbolic methods were never intended.
Performing numerical calculus requires a computer.
Less obviously, it also requires a way to communicate
with the computer, to tell the computer what to do. In
other words, to do calculus numerically, you need a way
to translate between human intention and electronic com-
putation: a language or notation.
These notes introduce one such notation for numerical
calculus, based on the computer language R.
R is a language for communicating instructions to a

computer and, it turns out, is also effective for communi-
cating with people. The R software is a computer system
that understands this language and acts on it.
Many people get nervous when they hear they will be
learning a new language. Most of our interaction with
computers word processing, e-mail, blogs and social
networking, spreadsheets is done with software thats
menu- and mouse-driven. Typically, its pretty easy to
learn such systems. You just have to be shown how and it
takes a few minutes to get started.
You will be spending a few hours getting started in
R. Why? Isnt there some simpler way? Why isnt there
nice mouse-based software? Why do you need to learn a
language?
The answer isnt about the quality of software or the
availability of friendly packages. Instead, the reason to
learn a language has to do with the sorts of things you
will be doing in mathematics and statistics. In them-
selves, the individual tasks you will undertake are not
necessarily more complicated than, say, changing a word
to a bold-face font, something you would do easily with a
mouse.
6
There are several aspects of technical comput-

ing that makes it different from word-processing and
other familiar, everyday computing tasks.
1. In mathematics and statistics, there are often multiple
inputs to a computation. To illustrate multiple inputs,
consider a familiar word-processing computation: find-
ing all the instances of the word car in a document
and changing them to automobile. Easy enough; just
use the find feature. But be careful! You might end
up with automobileeful or inautomobilecerated or
automobilecinogenic instead of careful, incarcerated,
or carcinogenic. A second input to the calculation is
needed: the set of contexts in which to allow or disal-
low the change. For instance, allow the change when
car is preceded by a space and followed by a space
or a period or a comma, or an s + space (cars !
automobiles). Things are not so simple as they might
seem at first, which is why the find-and-replace feature
of word-processors is only partially effective.
2. In mathematics and statistics, the output of one com-
putation often becomes the input to another compu-
tation. Thats why math courses spend so much time
talking about functions (and domain and range,
etc.). In word processing, whenever you highlight a
word and move it or change the font or replace it, you
still end up with stuff on which you can perform the
same operations: highlighting, moving, font-changing,
etc. Not so in math and statistics. The sorts of opera-
tions that you will often perform solving, integra-
tion, statistical summaries, etc. produce a new kind
of thing on which you will be performing new kinds
of operations. In mathematics and statistics, you cre-
ate a chain of operations and you need to be able to
express the steps in that chain. Its not a question of
having enough buttons to list all the operations, youll
need combinations of operations more than could
possibly be listed in a menu system.
3. In mathematics and statistics, the end-product is not
the only thing of importance. When you write a let-
ter or post to a blog, what counts is the final product,
not the changes you made while writing and certainly
7
not the thought process that you went through in com-

posing your words. But in mathematics and statistics,
the end-product is the result of a chain of calculations
and its important that each step in that chain be cor-
rect. Therefore, its important that each step in the
chain be documented and reproducible so that it can
be checked, updated, and verified. Often, the chain
of calculations becomes a new computation that you
might want to apply to a new set of inputs. Expressing
your calculations as a language allows you to do this.
If you have ever travelled to a country where you dont

speak the language, you know that you can communicate
simple ideas with gestures and pointing and can satisfy
the relatively simply expressed needs of eating and hy-
giene and shelter. But when you want to converse with
a person and express rich ideas, you need a shared lan-
guage.
Many people think that it would be better if other peo-
ple learned our language, and the natural extension of
this is that computers should be taught to speak English.
But it turns out that English, or other natural languages,
are not set up to be effective at communicating mathe-
matical or statistical ideas. You need to learn a way to
do this. The algebraic notation taught in high school is
part of the story, but not a complete solution. Thats why
youll be learning R.
If youve ever learned just a little of a foreign language,
youre familiar with the situation where you say some-
thing that seems straightforward, but your listener gives
you a quizzical look: it doesnt make sense. You used the
wrong verb or the wrong preposition or a word with a
slightly different meaning. A relative of mine, visiting me
in France, once asked my host about last year, or so he
intended. He actually said something pretty close to el-
der buttock. This did not produce the intended reaction.
Similarly in R. At first, you will make elementary mis-
takes. The computer will respond, like my French host,
quizzically. But with practice just a few hours you
will become fluent and able to express your ideas with
confidence and certainty.
As you learn the R language, which will be much,
much easier than learning a natural language like French
8
or Chinese or Spanish, you will make mistakes and you

will run into frustrating situations. But remember, the
reason you are learning it is to be able to express compli-
cated ideas. Its the nature of mathematics and statistics
thats at the core here. Having a systematic way to ex-
press yourself will not only let you use the computers
power, but will increase your understanding of the math-
ematical and statistical ideas.
My thanks go to my colleagues at Macalester, par-

ticularly Dan Flath and Chad Topaz, and the scores of
students from Math 135 who tested this material and
provided many suggestions for improvements. Randall
Pruim and Nicholas Horton have been instrumental in
transforming sketchy software ideas into the mosaic pack-
age. This work was supported in part by Project MOSAIC
and, especially, a grant from the US National Science
Foundation (DUE-0920350).
Daniel Kaplan, Dec. 2012
Contents
RStudio 11
Functions & Graphing 16
Fitting 34
Solving 41
Derivatives 54
Anti-Derivatives 59
Dynamics 64
Activities 69
Instructor Notes 73
10
Index 80
1. Starting with RStudio
R is a language that you can use to direct a computer to

perform mathematical and statistical operations, produce
graphics, carry out data processing. In between the
language and the computer is a system that interprets the
language and communicates your instructions to the
computers operating system and hardware. Due to the
popularity of R, there are many such interpreters. The
one you will be using is called RStudio.
You will likely use RStudio as a web service, like
Facebook or Google Docs.1 You need only have an 1
Its also possible to install
Internet connection and a recent web browser the sort RStudio on your own computer,
of thing you might use for Facebook. No other software and to use it without the web.
installation is required.
Figure 2: The RStudio window

displayed in a web browser.
The RStudio window has the familiar menu bar and is

divided into four panes:
1. The Console pane, into which you will type your
12 start r in calculus
commands. The console maintains a transcript of your

session with RStudio your command inputs and the
computers response.
2. The Workspace/History/Plots pane, which contains
three tabs. The Plots tab is where your graphics will
appear. History maintains an organized record of your
previous commands to help you remember.
3. The Files/Packages/Help pane, which allows you to
access documentation of commands as well as to load
in new specialized software (packages).
4. The Source pane, barely visible in the figure, which
you will not be using at first. Among other things,
source provides an editor for writing computer
programs, which are chains of commands stored for
later re-use.
Ready? Go to the console and give your first
command, right after the prompt >.
3+2
When you press return, R interprets your command

and gives its response.
[1] 5
Simple arithmetic in R is done with a familiar

notation. (See Figure 3 for some examples.) Remember to
use * multiplication and ^ for exponentiation.
Multiplication must always be specified explicitly. Its
correct to say 6*(3+1) but invalid to say 6(3+1).
rstudio 13
5-2 3^2 9*6+1
[1] 3 [1] 9 [1] 55
5/2 9*7 9*(6+1)
[1] 2.5 [1] 63 [1] 63
2^3
[1] 8
Figure 3: Arithmetic in R
For functions with letter names, for instance cos and

ln, you put the input to the function between parentheses.
For instance:
sqrt(9) acos(0) log(10)
[1] 3 [1] 1.571 [1] 2.303
cos(1) sin(pi) exp(-3)
[1] 0.5403 [1] 1.225e-16 [1] 0.04979
Figure 4: Use parentheses to

provide the input to named
functions.
Often, you will want to store the result of a
calculation under a name so that you can re-use the result
later. Such storage is called assignment and is
accomplished with the = operator:
first = 3
second = 4
hypot = sqrt(first^2 + second^2)
myangle = asin(second/hypot)
hypot*cos(myangle)
[1] 3
When you want to see the value of a named value, just

type the name as a command:
hypot
[1] 5
Of course, R can do much more than such

calculations, but thats enough to start for now.
Important Note about Typing: In the RStudio

console, you can use the keyboard up- and down-arrows
to recall previous statements. That saves you some typing
when you want to repeat something, or to make a small
change to a previous command. Remember that the
console is in the form of a dialog, not a document. Each
new command, even if it is an edited version of a
previous command, appears on the last line. The up- and
down-arrows copy previous commands to the current
prompt line, but you cant edit the past, just copy it.
rstudio 15
Exercises
Exercise 1 Use R to calculate a numerical Exercise 3 Here is a set of R expressions

value for each of these arithmetic that you should cut and paste as a command
expressions: into your R session:
alpha = 7
(a) 3*sqrt(97) beta = 3
13.72 21.61 23.14 29.55 31.01 31.20 33.03 x = 9
From these named quantities, you are
(b) sqrt(pi) going to compute a new one:
1.6152 1.7693 1.7725 1.8122 3.1416 3.519 res=atan(sqrt(alpha^2-beta)+
tan(x/alpha+beta/x))
(c) exp(pi)
13.72 21.61 23.14 29.55 31.01 31.20 33.03 You dont have to worry about what the
last expression might possibly mean, its just
intended to be something complicated. Make
(d) pi^3
sure to cut and paste it into R so that you get
13.72 21.61 23.14 29.55 31.01 31.20 33.03
it exactly right.
The value of res will be 1.4991. Before
Exercise 2 Use R to calculate a numerical continuing, make sure that your value
value for each of these arithmetic matches this, otherwise you wont get correct
expressions: answers for the following questions.
(a) Change alpha to 17 and recompute res.
p Whats the new value of res? (Hint: Use
(a) 2/3
the up arrow to recall the complicated
0.4805 0.5173 0.5642 0.6304 0.7071 0.8165
command.)
1.5099 1.5127 1.5155 1.5177 1.5182 1.5207
(b) cos p/4
0.4805 0.5173 0.5642 0.6304 0.7071 0.8165 (b) Keep alpha at 17 and x at 9, but change
beta to 6.6. Recompute res. Whats the
p new value of res?
(c) 1/ 2
0.4805 0.5173 0.5642 0.6304 0.7071 0.8165 1.5099 1.5127 1.5155 1.5177 1.5182 1.5207
(c) Keep alpha at 17 and beta at 6.6, but
p change x to 1. Recompute res. Whats the
(d) 1/ p
0.4805 0.5173 0.5642 0.6304 0.7071 0.8165 new value of res?
1.5099 1.5127 1.5155 1.5177 1.5182 1.5207
(e) 10 1/p (d) Keep alpha at 17 and beta at 6.6, but
0.4805 0.5173 0.5642 0.6304 0.7071 0.8165 change x to 1. Recompute res. Whats
the new value of res?
1.5099 1.5127 1.5155 1.5177 1.5182 1.5207
2. Functions & Graphing
2.1 Graphing Mathematical Functions

In this lesson, you will learn how to use R to graph
mathematical functions. The MOSAIC package.
Its important to point out at the beginning that much You will be using a hand-
of what you will be learning much of what will be new ful of functions in R, such as
plotFun() and makeFun(),
to you here actually has to do with the mathematical many of which are provided by
structure of functions and not R. the mosaic add-on package. To
Recall that a function is a transformation from an use these functions, you must
input to an output. Functions are used to represent the first tell R to load the package.
This can be done with this R
relationship between quantities. In evaluating a function,
command.
you specify what the input will be and the function
translates it into the output. require(mosaic)
In much of the traditional mathematics notation you
have used, functions have names like f or g or y, and the If you get an error message,
it is likely because the package
input is notated as x. Other letters are used to represent
has not yet been installed on
parameters. For instance, its common to write the your computer. Doing so is
equation of a line this way easy:
y = mx + b. install.packages("mosaic")
In order to apply mathematical concepts to realistic You need install the package
only once. (It may already
settings in the world, its important to recognize three have been done for you.) But
things that a notation like y = mx + b does not support you need to load the mosaic
well: package each time you restart
R.
1. Real-world relationships generally involve more than
two quantities. (For example, the Ideal Gas Law in
chemistry, PV = nRT, involves three variables:
pressure, volume, and temperature.) For this reason,
you will need a notation that lets you describe the
multiple inputs to a function and which lets you keep
track of which input is which.
functions & graphing 17
2. Real-world quantities are not typically named x and y,

but are quantities like cyclic AMP concentration or
membrane voltage or government expenditures.
Of course, you could call all such things x or y, but its
much easier to make sense of things when the names
remind you of the quantity being represented.
3. Real-world situations involve many different
relationships, and mathematical models of them can
involve different approximations and representations
of those relationships. Therefore, its important to be
able to give names to relationships, so that you can
keep track of the various things you are working with.
For these reasons, the notation that you will use needs
to be more general than the notation commonly used in
high-school algebra. At first, this will seem odd, but the
oddness doesnt have to do so much with the fact that the
notation is used by the computer so much as for the
mathematical reasons given above.
But there is one aspect of the notation that stems
directly from the use of the keyboard to communicate
with the computer. In writing mathematical operations,
youll use expressions like a*b and 2^n and a/b rather
than the traditional ab or 2n or ba , and you will use
parentheses both for grouping expressions and for
applying functions to their inputs.
In plotting a function, you need to specify several
things:
What is the function. This is usually given by an
expression, for instance m*x + b or A*x^2 or sin(2*t)
Later on, you will also give names to functions and use
those names in the expressions, much like sin is the
name of a trigonometric function.
What are the inputs. Remember, theres no reason to
assume that x is always the input, and youll be using
variables with names like G and cAMP. So you have to
be explicit in saying whats an input and whats not.
The R notation for this involves the ~ (tilde) symbol.
For instance, to specify a linear function with x as the
input, you can write m*x + b ~ x
What range of inputs to make the plot over. Think of this as
the bounds of the horizontal axis over which you want

to make the plot.
The values of any parameters. Remember, the notation
m*x+b ~ x involves not just the variable input x but
also two other quantities, m and b. To make a plot of
the function, you need to pick specific values for m and
b and tell the computer what these are.
The plotFun() operator puts this all together, taking 25

the information you give and turning it into a plot. Heres 20
3*x-2
an example of plotting out a linear function: 15
10
plotFun(3*x - 2 ~ x, x.lim=range(0,10) ) 5
Often, its natural to write such relationships with the 2 4 6 8

parameters represented by symbols. (This can help you x
remember which parameter is which, e.g., which is the
slope and which is the intercept. When you do this,
remember to give a specific numerical value for the
parameters, like this:
25
20
plotFun(m*x + b ~ x, x.lim=range(0,10), m=3, b=-2 )
m*x+b
15
10
Try these examples: 5
0
plotFun(A*x^2 ~ x, x.lim=range(-2,3), A=10)
plotFun(A*x^2 ~ x, add=TRUE, col="red", A=5) 2 4 6 8
x
plotFun(cos(t) ~ t, t.lim=range(0,4*pi))
Sometimes, you want to give a function a name so that

you can refer to it concisely later on. You can use
makeFun() to create a function and ordinary assignment
to give the function a name. For instance:
g = makeFun(2*x^2 - 5*x + 2 ~ x)
Once the function is named, you can evaluate it by

giving an input. For instance:
g(x=2) g(x=5)
[1] 0 [1] 27
To plot g(), you need to provide an argument and use

the ~ to identify which variable corresponds to the x-axis.
plotFun(g(x) ~ x, x.lim=range(-5,5))
Exercises
Exercise 1 This command Exercise 4

plotFun(A*x^2 ~ A, A.lim=range(-2,3), x=10) Reproduce each of these two plots. Hand
in the command you used to make the
will graph a straight line. Explain why identical plots:
the graph doesnt look like a parabola, even
though its a graph of Ax2 . 4
6
3
Exercise 2 4
2*x3
t^2
2 2
Translate each of these expressions in 0
traditional math notation into a plot made by
1
2
plotFun(). Hand in the command that you

0
1 2 3 4 1 0 1
gave to make the plot (not the plot itself). x t
(a) 4x 7 in the window x from 0 to 10.

(b) cos 5x in the window x from 1 to 1.
(c) cos 2t in the window t from 0 to 5. Exercise 5
p What happens when you use a symbolic
(d) t cos p 5t in the window t from 0 to 5.
parameter (e.g., m in m*x + b ~ x, but try to
(Hint: (t) is sqrt(t).)
make a plot without selecting a specific
numerical value for the parameter?
Exercise 3 Find the value of each of the
functions in Exercise 2 at x = 10.543 or at
t = 10.543. (Hint: Give the function a name Exercise 6
and compute the value using an expression What happens when you dont specify a
like g(x=10.543) or f(t=10.543).) range for an input, but just a single number,
Pick the closest numerical value as in the second of these two commands:
(a) 32.721 34.721 35.172 37.421 37.721 plotFun(3*x ~ x, x.lim=range(1,4))
(b) -0.83 -0.77 -0.72 -0.68 0.32 0.42 0.62 plotFun(3*x ~ x, x.lim=14)
(c) -0.83 -0.77 -0.72 -0.68 -0.62 0.42 0.62 Give a description of what happened and
(d) -2.5 -1.5 -0.5 0.5 1.5 2.5 speculate on why.
2.2 Making Scatterplots from Data

Often, the mathematical models that you will create will
be motivated by data. For a deep appreciation of the
relationship between data and models, you will want to
study statistical modeling. Here, though, we will take a
first cut at the subject in the form of curve fitting, the
process of setting parameters of a mathematical function
to make the function a close representation of some data.
This means that you will have to learn something
about how to access data in computer files, how data are
stored, and how to visualize the data. Fortunately, R and
the mosaic package make this straightforward.
The data files you will be using are stored as
spreadsheets on the Internet. Typically, the spreadsheet
will have multiple variables; each variable is stored as one
column. (The rows are cases, sometimes called data
points.) To read the data in to R, you need to know the
name of the file and its location.
The function for reading such files from their
appointed place on the Internet is called fetchData().
One of the files that fetchData() can read is
"Income-Housing.csv". This file gives information from a
survey on housing conditions for people in different
income brackets in the US.2 Heres how to read it into R: 2
Susan E. Mayer (1997) What
money cant buy: Family income
housing = fetchData("Income-Housing.csv") and childrens life chances Har-
vard Univ. Press p. 102.
There are two important things to notice about the
above statement. First, the fetchData() function is
returning a value that is being stored in an object called
housing. The choice of housing as a name is arbitrary;
you could have stored it as x or Equador or whatever. Its
convenient to pick names that help you remember whats
being stored where.
Second, the name "Income-Housing.csv" is
surrounded by quotation marks. These are the
single-character double quotes, that is, " and not repeated
single quotes . Whenever you are reading data from a
file, the name of the file should be in such
single-character double quotes. That way, R knows to
treat the characters literally and not as the name of an
object such as housing.
Once the data are read in, you can look at the data just
by typing the name of the object (without quotes!) that is
holding the data. For instance,
housing
Income IncomePercentile CrimeProblem AbandonedBuildings

1 3914 5 39.6 12.6
2 10817 15 32.4 10.0
3 21097 30 26.7 7.1
4 34548 50 23.9 4.1
5 51941 70 21.4 2.3
6 72079 90 19.9 1.2
All of the variables in the dataset will be shown

(although just four of them are printed here).
You can see the names of all of the variables in a
compact format with the names() command:
names(housing)
[1] "Income" "IncomePercentile" "CrimeProblem"

[4] "AbandonedBuildings" "IncompleteBathroom" "NoCentralHeat"
[7] "ExposedWires" "AirConditioning" "TwoBathrooms"
[10] "MotorVehicle" "TwoVehicles" "ClothesWasher"
[13] "ClothesDryer" "Dishwasher" "Telephone"
[16] "DoctorVisitsUnder7" "DoctorVisits7To18" "NoDoctorVisitUnder7"
[19] "NoDoctorVisit7To18"
When you want to access one of the variables, you

give the name of the whole dataset followed by the name
of the variable, with the two names separated by a $ sign,
like this:
housing$Income
[1] 3914 10817 21097 34548 51941 72079
housing$CrimeProblem
[1] 39.6 32.4 26.7 23.9 21.4 19.9
Even though the output from names() shows the

variable names in quotation marks, you wont use
quotations around the variable names.
Spelling and capitalization are important. If you make

a mistake, no matter how trifling to a human reader, R
will not figure out what you want. For instance, heres a
misspelling of a variable name, which results in nothing
(that is, NULL) being returned.
housing$crim
NULL
Sometimes people like to look at datasets in a

spreadsheet format, each entry in a little cell. In RStudio,
you can do this by going to the Workspace tab and
clicking the name of the variable you want to look at.
Figure 5: Viewing a dataset by

clicking on the object name in
the Workspace tab in RStudio.
40
CrimeProblem
35
Its hard to see patterns in columns of numbers, but a 30
graphical presentation of data can be informative. One of 25

the most familiar graphical forms is the scatterplot, a
format in which each "case" or "data point" is plotted as a 20
dot at the coordinate location given by two variables. For 20000 40000 60000
Income
instance, heres a scatter plot of the percentage of
household that regard their neighborhood as having a Figure 6: A scatterplot.
crime problem, versus the median income in their
bracket.
plotPoints( CrimeProblem ~ Income, data=housing ) 40

CrimeProblem
35
The R statement closely follows the English
equivalent: plot CrimeProblem versus (or, as a function 30
of) Income, using the data from the housing object. 25

If you want to plot a mathematical function over the
data, do so by first giving a plot command to show the
20
data, then asking plotFun() to add a graph of the 20000 40000 60000
Income
function:
Figure 7: Adding a mathemati-
cal function to the scatterplot.
plotPoints( CrimeProblem ~ Income, data=housing )

plotFun( 40 - Income/2000 ~ Income,
Income.lim=range(0,80000), add=TRUE )
The function drawn is not a very good match to the

data, but this reading is about how to draw graphs, not
how to choose a family of functions or find parameters!
Notice the add=TRUE argument to plotFun(), which
instructs R to add the new graph over the old one.
Without this, R would draw a brand-new graph.
The plotFun() graph-drawing function allows you to
give your mathematical function arguments of whatever
name you like. So you could add another graph to the
plot by giving a function like this:
plotFun( 38 - x/3500 ~ x, x.lim=range(0,80000), add=TRUE, col="red")
If, when plotting your data, you prefer to set the limits
of the axes to something of your own choice, you can do
this. For instance:
40
plotPoints( CrimeProblem ~ Income, data=housing,
CrimeProblem
xlim=range(0,100000), ylim=range(0,50) ) 30
plotFun( 40 - Income/2000 ~ Income, 20
Income=range(0,80000), add=TRUE ) 10
Properly made scientific graphics should have 2e+04 4e+04 6e+04 8e+04
informative axis names. You can set the axis names Income
directly in either plotFun() or plotPoints():
plotPoints( CrimeProblem ~ Income/1000, data=housing, Crime Problem

xlab="Income Bracket (USD/year)",
ylab="Percentage of Households",
Percentage of Households
40
main="Crime Problem",
30
xlim=range(0,100), ylim=range(0,45) )
20
Notice that double-quotes delimit the character 10

strings. The argument names xlim, ylim, xlab and ylab
are used to refer to the ranges and labels of the horizontal 20 40 60 80
and vertical axes respectively. Income Bracket (USD/year)

Exercises
Exercise 1 Exercise 2
Make each of these plots: Construct the R commands to duplicate
(a) Prof. Stan Wagon (see each of these plots. Hand in your commands
http://stanwagon.com) illustrates curve (not the plot):
fitting using measurements of the (a) The data file "utilities.csv" has utility
temperature (in degrees C) of a cup of records for a house in St. Paul, Minnesota,
coffee versus time (in minutes): USA. Make this plot, including the labels:
Ave. Monthly Temp.
s = fetchData("stan-data.csv")
plotPoints(temp ~ time, data=s) 80
Temperature (F)
60
100 40
80 20
temp
60 2 4 6 8 10 12
Month (Jan=1, Dec=12)

40
0 50 100 150 200

(b) From the "utilities.csv" data file, make
time this plot of household monthly bill for
natural gas versus average temperature.
Describe in everyday English the pattern
The line has slope 5 $US/degree and
you see in coffee cooling:
intercept 300 $US.
(b) Heres a record of the tide level in Hawaii Natural Gas Use
over about 100 hours:
Expenditures (USD per month)
250
h = fetchData("hawaii.csv") 200
plotPoints(water ~ time, data=h) 150
100
50
0
2.0
20 40 60 80
1.5 Temperature (F)
water
1.0
0.5
0.0
0 20 40 60 80
time
Describe in everyday English the pattern

you see in the tide data:
2.3 Graphing Functions of Two Variables

Youve already seen how to plot a graph of a function of
95 - 73 * exp(-0.2 * t)
one variable, for instance: 80
60
plotFun(95-73*exp(-.2*t) ~ t,
t.lim=range(0,20)) 40
20
This lesson is about plotting functions of two 5 10 15
variables. For the most part, the format used will be a t
contour plot, but its also possible to make the graph of
the function, as youll see later.
You use the same plotFun() function to plot with two
input variables. The only change is that you need to list
the two variables on the right of the ~ sign, and you need
to give a range for each of the variables. For example:
0.0
8
plotFun( sin(2*pi*t/10)*exp(-.2*x) ~ t & x,
t.lim=range(0,20),x.lim=range(0,10)) 6
x
4 -0.5 0.5
0.5 -0.5
Each of the contours is labeled, and by default the plot
0.0
0.0
2
is filled with color to help guide the eye. If you prefer just
to see the contours, without the color fill, use the 5 10 15
filled=FALSE argument. t
plotFun( sin(2*pi*t/10)*exp(-.2*x) ~ t & x,

t.lim=range(0,20),x.lim=range(0,10),filled=FALSE)
0.0
Occasionally, people want to see the function as a 8
surface, plotted in 3 dimensions. You can get the
6
computer to display a perspective 3-dimensional plot by
x
giving the optional argument surface=TRUE. 4

0.5 -0.5 0.5 -0.5
0.0
0.0
plotFun(sin(2*pi*t/10)*exp(-.2*x) ~ t & x,
t.lim=range(0,20),x.lim=range(0,10),surface=TRUE) 5 10 15
t
If you are using RStudio, you can press on the little
gear icon in the plot and you will have a slider to
control the viewpoint. (Try moving the slider to the right,
release it, and wait for the picture to update.)
Its very hard to read quantitative values from a
surface plot the contour plots are much more useful 0.5
for that. On the other hand, people seem to have a strong 0.0
-0.5
10
8 20
6 15
4 10
x 2 5 t
00
intuition about shapes of surfaces. Being able to translate

in your mind from contours to surfaces (and vice versa) is
a valuable skill.
Sometimes you will want to create and plot a

named function. When the function has multiple
arguments, you need to be aware of which argument is
which when you plot it, or you might get them reversed.
Heres a safe strategy. 0.0
8
g = makeFun( sin(2*pi*t/10)*exp(-.2*x) ~ t & x) 6
x
plotFun(g(t=t,x=x) ~ t & x, 4
0.5 -0.5 0.5 -0.5
0.0
0.0
t.lim=range(0,20),x.lim=range(0,10)) 2
The point of the seemingly redundant t=t and x=x is 5 10 15

to avoid mistaking the arguments one for the other. The t
name to the left of the = sign gives the name of the input
variable. The expression g(t=t,x=x) says to assign the
plotting variable t (which will range from 0 to 20) to the
input named t.
So long as you use the input names, the order doesnt
matter:
g(t=7,x=4) g(x=4,t=7)
[1] -0.4273 [1] -0.4273
But if you dont use the names, you might end up

assigning the input value to the wrong input:
g(7,4) g(4,7)
[1] -0.4273 [1] 0.1449

Exercises
Refer to this contour plot: Describe the shape of the contours
produced by each of these functions. (Hint:
6 Make the plot! Caution: Use the mouse to
6
-6
4
make the plotting frame more-or-less square
4
2
8
2
-4
6
in shape.)
x
4 -2 0
(a) The function
0
plotFun( sqrt((v-3)^2+2*(w-4)^2) ~ v & w,

5 10 15 v.lim=range(0,6),w.lim=range(0,6))
t
has contours that are
Approximately what is the value of the A Parallel Lines
function at each of these ( x, t) pairs? Pick the B Concentric Circles
closest value C Concentric Ellipses
D X Shaped
(b) The function
(a) x = 4, t = 10: -6 -5 -4 -2 0 2 4 5 6 plotFun( sqrt( (v-3)^2+(w-4)^2 ) ~ v & w,

v.lim=range(0,6),w.lim=range(0,6))
has contours that are
A Parallel Lines
(b) x = 8, t = 10: -6 -5 -4 -2 0 2 4 5 6 B Concentric Circles
C Concentric Ellipses
D X Shaped
(c) The function
(c) x = 7, t = 0: -6 -5 -4 -2 0 2 4 5 6
plotFun( 6*v-3*w+4 ~ v & w,
v.lim=range(0,6),w.lim=range(0,6))
has contours that are:
(d) x = 9, t = 0: -6 -5 -4 -2 0 2 4 5 6 A Parallel Lines
B Concentric Circles
C Concentric Ellipses
D X Shaped
2.4 Splines & Smoothers

Mathematical models attempt to capture patterns in the
real world. This is useful because the models can be more
easily studied and manipulated than the world itself. One
of the most important uses of functions is to reproduce or
capture or model the patterns that appear in data.
Sometimes, the choice of a particular form of function
exponential or power-law, say is motivated by an
understanding of the processes involved in the pattern
that the function is being used to model. But other times,
all thats called for is a function that follows the data and
that has other desirable properties, for example is smooth
or increases steadily.
Smoothers and splines are two kinds of
general-purpose functions that can capture patterns in
data, but for which there is no simple algebraic form.
Creating such functions is remarkably easy, so long as
you can free yourself from the idea that functions must
always have algebraic formulas like f ( x ) = ax2 + b .
Smoothers and splines are defined not by algebraic
forms and parameters, but by data and algorithms. To
illustrate, consider some simple data. The data set
Loblolly contains 84 measurements of the age and height
of loblolly pines.
60
pine = fetchData("Loblolly") 50
plotPoints(height ~ age,data=pine)
height
40
30
Several three-year old pines of very similar height 20
were measured and tracked over time: age five, age ten, 10
and so on. The trees differ from one another, but they are
5 10 15 20 25
all pretty similar and show a simple pattern: linear
age
growth at first which seems to slow down over time.
It might be interesting to speculate about what sort of
algebraic function the loblolly pines growth follows, but
any such function is just a model. For many purposes,
measuring how the growth rate changes as the trees age,
all thats needed is a smooth function that looks like the
data. Lets consider two:
A cubic spline, which follows the groups of data

points and curves smoothly and gracefully.
f1 = spliner(height ~ age, data=pine)
A linear interpolant, which connects the groups of

data points with straight lines.
f2 = connector(height ~ age, data=pine)
The definitions of these functions may seem strange at

first they are entirely defined by the data: no
parameters! Nonetheless, they are genuine functions and
can be worked with like other functions. For example,
you can put in an input and get an output:
f1(age=8)
[1] 20.68
f2(age=8)
[1] 20.55
You can graph them: 60
plotFun(f1(age) ~ age, age.lim=range(0,30)) f1(age) 40
plotFun(f2(age) ~ age, add=TRUE, col="red")

20
plotPoints(height ~ age, data=pine, add=TRUE)
0
In all respects, these are perfectly ordinary functions. 5 10 15 20 25

All respects but one: There is no simple formula for them. age
Youll notice this if you ever try to look at the
computer-language definition of the functions:
f2
function (age)
{
x <- get(fnames[2])
if (connect)
SF(x)
else SF(x, deriv = deriv)
}
<environment: 0x1046c06d8>
Theres almost nothing here to tell you what the

function is. The definition refers to the data itself which
has been stored in an environment. These are
computer-age functions, not functions from the age of
algebra.
As you can see, the spline and linear connector
functions are quite similar, except for the range of inputs
outside of the range of the data. Within the range of the
data, however, both types of functions go exactly through
the center of each age-group.
Splines and connectors are not always what you will
want, especially when the data are not divided into
discrete groups, as with the loblolly pine data. For
instance, the trees data set is measurements of the
volume, girth, and height of black cherry trees. The trees
were felled for their wood, and the interest in making the
measurements was to help estimate how much usable
volume of wood can be gotten from a tree, based on the
girth (that is, circumference) and height. This would be
useful, for instance, in estimating how much money a tree
is worth. However, unlike the loblolly pine data, the black
cherry data does not involve trees falling nicely into
80
defined groups. Volume 60
cherry = fetchData("trees") 40
plotPoints(Volume~Girth, data=cherry)
20
Its easy enough to make a spline or a linear 10 15 20

connector: Girth
g1 = spliner(Volume~Girth, data=cherry)
g2 = connector(Volume~Girth, data=cherry)
plotFun(g1(x)~x,x.lim=range(8,18),xlab="Girth (inches)")
plotFun(g2(x) ~ x, add=TRUE, col="red")
plotPoints(Volume ~ Girth, data=cherry, add=TRUE)
The two functions both follow the data ... but a bit too
faithfully! Each of the functions insists on going through 60
every data point. (The one exception is the two points 50
with girth of 13 inches. Theres no function that can go 40

g1(x)
through both of the points with girth 13, so the functions 30
split the difference and go through the average of the two 20

points.) 10
10 12 14 16
Girth (inches)
Its hard to believe that the rapid up-and-down

wiggling is of the functions is realistic. When you have
reason to believe that a smooth function is more
appropriate than one with lots of ups-and-downs, a 50
different type of function is appropriate: a smoother. 40
g3(x)
30
g3 = smoother(Volume ~ Girth, data=cherry, span=0.5) 20
plotFun(g3(x)~x,x.lim=range(8,18),xlab="Girth (inches)") 10
plotPoints(Volume~Girth, data=cherry, add=TRUE) 10 12 14 16
Girth (inches)
Smoothers are well named: they construct a smooth
function that goes close to the data. You have some
control over how smooth the function should be. The
parameter span governs this: 50
40
g4(x)
g4 = smoother(Volume ~ Girth, data=cherry, span=1.0) 30
20
Of course, often you will want to capture relationships 10
where there is more than one variable as the input. 10 12 14 16
Smoothers do this very nicely; just specify which Girth (inches)
variables are to be the inputs.
Figure 8: Control the smooth-
ness of the curve with the span
g5 = smoother(Volume ~ Girth+Height, data=cherry,
parameter.
span=1.0)
plotFun(g5(g,h) ~ g & h,
g.lim=range(8,18), h.lim=range(60,90))
When you make a smoother or a spline or a linear

connector, remember these rules:
You need a dataset containing the data.

50
85
80
You use the formula with the variable you want as the 75
h
40
output of the function on the left side of the tilde, and 70
20
30
the input variables on the right side. 65

10
The function that is created will have input names that 10 12 14 16

match the variables you specified as inputs. (For the g
present, only smoother will accept more than one input
Figure 9: A smoother in two
variable.)
variables: Girth and Height.
The smoothness of a smoother function can be set by
the span argument. A span of 1.0 is typically pretty
smooth. The default is 0.5.
When creating a spline, you have the option of

declaring monotonic=TRUE. This will arrange things to
avoid extraneous bumps in data that shows a steady
upward pattern or a steady downward pattern.
When you want to plot out a function, you need of

course to choose a range for the input values. Its often
sensible to select a range that corresponds to the data on
which the function is based. You can find this with the
max() and min() commands, e.g.
min(Height, data=cherry) max(Height, data=cherry)
[1] 63 [1] 87
Exercises
Exercise 1 A The span=10 smoother is

The Orange data set contains data on the always above the span=1
circumference and age of trees. smoother.
orange = fetchData("Orange") B The span=10 smoother is
names(orange)
always below the span=1
smoother.
[1] "Tree" "age"
C The span=10 smoother is
[3] "circumference"
straighter than the span=1
smoother.
(a) Construct a smoother of circumference (in
cm) as a function of age (in days), with
span=1. Whats the value of the function Exercise 2
for age 1000 days? The cherry-tree data gives wood volume
127.41 128.32 130.28 132.08 133.12 138.20 (in cubic feet) along with height (in feet) and
girth (in inches even though the name
(b) Construct a smoother of circumference Girth is used, its really a diameter).
versus age with span=10. Whats the Construct a smoother of two variables for
value for age 1000 days? the cherry data with span=1. (Warning: Be
127.41 128.32 130.28 132.08 133.12 138.20 careful to match the capitalization of the
names).
(c) How does the shape of the graph of the Whats the value when height is 50 feet
span=10 smoother differ from that of the and Girth (that is, diameter) is 15 inches?
span=1 smoother? 24.4 25.5 26.6 27.7 28.8
Exercise 3 23.34 23.85 24.43 25.23 26.74 27.25

The data set BodyFat.csv contains
(ii) Whats the general shape of the
several body measurements made on a
function?
group of 252 men. Youre going to look at A Sloping up and curved down-
theses variables: wards
Weight the total body weight in B Sloping down and curved
pounds. upwards
BodyFat the percentage of total weight C Straight line sloping down-
that is fat. wards
Abdomen the circumference at the waist
in cm. (c) Construct a smoother of BodyFat as a
function of both Abdomen and Weight with
(a) Construct a smoother of BodyFat as a
span=1. What is the value of this function
function of Weight with span=1.
for Abdomen=100 and Weight=150?
(i) What is the value of this function for 23.34 23.85 24.43 25.23 26.74 27.25
Weight=150? 13.46 14.36 16.34
(ii) Whats the general shape of the Make a contour plot of your smoother
function? function. Using the plot, answer this
A Sloping up and curved down- question: How does body fat percentage
wards change when increasing weight from 140
B Sloping down and curved to 240 pounds, but holding abdominal
upwards circumference constant at 90cm?
C Straight line sloping down- A Body fat percentage doesnt
wards change on weight.
(b) Construct a smoother of BodyFat as a B Body fat percentage goes down
function of Abdomen with span=1. with increasing weight.
(i) What is the value of this function for C Body fat percentage goes up
Abdomen=100? with increasing weight.
3. Fitting Functions to Data
Often, you have an idea for the form of a function for a

model and you need to select parameters that will make
the model function a good match for observations. The
process of selecting parameters to match observations is
called model fitting.
To illustrate, the data in the file "utilities.csv"
records the average temperature each month (in degrees
F) as well as the monthly natural gas usage (in cubic feet,
ccf). There is, as you might expect, a strong relationship
between the two.
250
u = fetchData("utilities.csv") 200
plotPoints(ccf ~ temp, data=u) 150

ccf
100
Many different sorts of functions might be used to 50
represent these data. One of the simplest and most
0
commonly used in modeling is a straight-line function
20 40 60 80
f ( x ) = Ax + B. In function f ( x ), the variable x stands for
temp
the input, while A and B are parameters. Its important to
remember what are the names of the inputs and outputs
when fitting models to data you need to arrange for
the name to match the corresponding data.
With the utilities data, the input is the temperature,
temp. The output that is to be modeled is ccf. To fit the
model function to the data, you write down the formula
with the appropriate names of inputs, parameters, and
the output in the right places:
f = fitModel(ccf ~ A*temp + B, data=u)
The output of fitModel() is a function of the same

form as you specified with specific numerical values
given to the parameters in order to make the function
best match the data.
fitting 35
The point of the fitting process is to find the

coefficients that bring the function as close as possible to
the data. You can easily see the coefficients that are found
by fitting by applying the coef() operator to the model.
coef(f)
A B
-3.464 253.098
Most often, you will be interested in the function

itself, for example to graph it along with the data.
250
plotFun( f(temp)~temp, temp.lim=range(0,80)) 200
plotPoints(ccf ~ temp, data=u, add=TRUE)
f(temp)
150
100
You can add other functions into the mix easily. For 50
instance, you might think that sqrt(temp) works in there 0
somehow. Try it out!

20 40 60
temp
f2 = fitModel(ccf ~ A*temp + B + C*sqrt(temp),
data=u)
plotFun( f2(temp)~temp, temp.lim = range(0,80))
plotPoints(ccf ~ temp, data=u, add=TRUE)
400
This example has involved just one input variable. f2(temp)

300
Throughout the natural and social sciences, a very 200
important and widely used technique is to use multiple
variables as inputs. To illustrate, look at the data in
100
"used-hondas.csv" on the prices of used Honda 0
automobiles. 20 40 60
temp
hondas = fetchData("used-hondas.csv")
head(hondas)
Price Year Mileage Location Color Age

1 20746 2006 18394 St.Paul Grey 1
2 19787 2007 8 St.Paul Black 0
3 17987 2005 39998 St.Paul Grey 2
4 17588 2004 35882 St.Paul Black 3
5 16987 2004 25306 St.Paul Grey 3
6 16987 2005 33399 St.Paul Black 2
As you can see, the dataset includes the variables

Price, Age, and Mileage. It seems reasonable to think that
price will depend both on the mileage and age of the car.
Heres a very simple model that uses both variables:
carPrice1 = fitModel(
Price ~ A + B*Age + C*Mileage,
data=hondas)
You can plot that out as a mathematical function:

13000
15000
plotFun( carPrice1(Age=age, Mileage=miles)~age&miles, 50000
14000
age.lim=range(2,8), miles.lim=range(0,60000)) 40000
miles
30000
18000
16000
A somewhat more sophisticated model might include 20000
17000
whats called an interaction between age and mileage, 10000
20 19000
00
recognizing that the effect of age might be different 0
3 4 5 6 7
depending on mileage.
age
carPrice2 = fitModel(Price ~ A+B*Age+C*Mileage+D*Age*Mileage,

data=hondas)
Again, once the function has been fitted, you can plot
it in the ordinary way:
10
00
0
plotFun( carPrice2(Age=age, Mileage=miles)~age & miles, 80000 12
00
0
age.lim=range(0,10), miles.lim=range(0,100000)) 60000
14
00
miles
16 0
40000 00
0
Notice that the price of a used car goes down with age 18
000
20
and with mileage. This is hardly unexpected. The fitted 20000 00
0
model quantifies the relationship, and from the graph you
2 4 6 8
can see that the effect of 2 years of age is roughly the
age
same as 20,000 miles.
Each of the above models has involved what are called
linear parameters. Often, there are parameters in
functions that appear in a nonlinear way. Examples
include k in f (t) = A exp(kt) + C and P in A sin( 2pP t ) + C.
The idea of function fitting applies perfectly well to
nonlinear parameters, but the task is harder for the
computer. Youll get the best results if you give the
computer a hint for the values of nonlinear parameters.
To illustrate, consider the "Income-Housing.csv" data
which shows an exponential relationship between the
fraction of families with two cars and income:
fitting 37
inc = fetchData("Income-Housing.csv")
plotPoints(TwoVehicles ~ Income, data=inc)
The pattern of the data suggests exponential decay 80
TwoVehicles
towards close to 100% of the families having two vehicles. 60
The mathematical form of this exponential function is

40
Aexp( kY ) + C. A and C are unknown linear
parameters. k is an unknown nonlinear parameter it 20
will be negative for exponential decay. 20000 40000 60000
Suppose you make a guess at k. The guess doesnt Income

need to be completely random; you can see from the data
themselves that the half-life is something like $25,000.
The parameter k is corresponds to the half life,
k = ln(0.5)/half-life, so a good guess for k is
ln(0.5)/25000, that is
kguess = log(0.5)/25000
kguess
[1] -2.773e-05
Its also helpful to have reasonable guesses for the

other parameters. You can make reasonable estimates
from the graph. A corresponds to the "leveling-off" value
of the function, which is around 100. B is negative and
about the same size as A (since essentially nobody has
two vehicles for a family income of 0). Starting with those
guesses, fitModel() can find good-fitting values for the
parameters:
f = fitModel( TwoVehicles ~ A + B*exp(k*Income), data=inc,

start=list(A=100,B=-100,k=log(0.5)/25000))
100
plotFun(f(Income)~Income, Income.lim=range(0,100000)) 80
plotPoints(TwoVehicles ~ Income, data=inc, add=TRUE)

f(Income)
60
40
The graph goes satisfyingly close to the data points.
20
But you can also look at the numerical values of the
function for any income: 0
2e+04 4e+04 6e+04 8e+04
Income
f(Income=10000)
[1] 33.78
f(Income=50000)
[1] 85.44
Its particularly informative to look at the values of the

function for the specific Income levels in the data used for
fitting, that is, the dataset inc:
f(Income=inc$Income)
[1] 16.38 35.82 56.67 74.01 86.45 93.48
The residuals are the difference between these model

values and the actual values of TwoVehicles in the
dataset:
resids = inc$TwoVehicles - f(Income=inc$Income)

resids
[1] 0.9247 -1.5186 -0.2713 1.2948 0.1478 -0.5774
This set of numbers is a vector; its length tells how

"far" the function is from the data. Recall that the
square-length of a vector is the sum of squared residuals
sum(resids^2)
[1] 5.267
The "distance" or rather the square distance

between the function and the data is the sum of square
residuals.
Keep in mind that the sum of square residuals is a 100
function of the parameters. The parameters are being 80

chosen by fitModel() in order to make the sum of square
temp
residuals as small as possible, in other words, to fit the 60
function to the data. 40
Example: The cooling of a hot object to the ambient 0 50 100 150 200
temperature is generally modeled by an exponential time
process. Lets see. The data "stan-data.csv" contain
Figure 10: The temperature
temperature measurements of a cup of hot water as it of cooling water and a fitted
cools to room temperature. To fit an exponential decay exponential model.
model, T = A + Be kt , well need an estimate for the
nonlinear parameter k. Looking at a plot of the data
fitting 39
suggests that it takes about 50 seconds for the initial

temperature to fall by about half of its eventual total fall.
This suggests k = ln(2)/50.
water = fetchData("stan-data.csv")
plotPoints(temp~time, data=water)
f = fitModel(temp ~ A + B*exp(-k*time), data=water,
start=list(k=log(2)/50))
plotFun( f(time)~time, add=TRUE, col="red")
You can see from the plot that the model captures the
gross shape of the data, but deviates from it at the start
and at the end. Its helpful to plot out the residuals the 8
data values minus their corresponding model values.
temp - f(time)
6
4
plotPoints( temp - f(time) ~ time, data=water) 2
Ideally, these should be near zero and show no -2
particular pattern in their deviations from zero. As you 0 50 100 150 200
can see, however, these residuals show systematic trends: time
slow oscillations above and below zero. That indicates
Figure 11: Residuals from the
that the model is not representing the cooling process
model show systematic trends
very well. It turns out that there are at least two cooling around zero.
processes, water to mug and water to air. A model with
two exponentials, one fast and one slow, does a much
better job.
Exercises
Exercise 1 The data in "stan-data.csv" (a) What is the period P (in hours) that
contains measurements made by Prof. Stan makes the sum of square residuals error
Wagon of the temperature of a cooling cup of as small as possible?
hot water. The variables are temp and time: 23.42 24.00 24.28 24.54 24.78 25.17
temperature in degrees C and time in
(b) Plot out the data and the fitted function.
minutes.
You may notice that the best fitting sine
Find the best value of k in the
wave is not particularly close to the data
exponential model A + B exp(kt).
points. One reason for this is that the
(a) Whats the value of k that gives the pattern is more complicated than a simple
smallest sum of square residuals? (Pick sine wave. You can get a better
the closest one.) approximation by including additional
-2.00 -0.20 -0.02 -0.002 -0.0002 sine functions with a period of 2P. (This
(b) What are the units of this k? (This is not is called a harmonic.) Overall, the model
an R question, but a mathematical one.) function will be:
A seconds
B minutes 2p
f 2(t) = A sin (t T0 ) +
C per second P

D per minute 2p
B sin (t T1 ) + C.
P/2
Exercise 2 The "hawaii.csv" dataset What period P (in hours) shows up as

contains a record of ocean tide levels in best when you add in a harmonic to the
Hawaii over a few days. The time variable is model?
in hours. You are going to fit the function 23.42 24.00 24.28 24.54 24.78 25.17
f (t) = A sin 2p
P (t T0 ) + C. Since the tides
occur with a period of roughly 1 day, a good A more complete model of tides includes
guess for a starting value of P is 24 hours. multiple periods stemming from the
multiple factors involved: the earths
hawaii = fetchData("hawaii.csv") rotation, the moons revolution around
f=fitModel(water~A*sin(2*pi*(time-T0)/P)+C, the earth, the alignment of the moon and
start=list(P=24,T0=0), data=hawaii) the sun, etc.
4. Solving
4.1 Solving Equations

Many of high-school algebra involves solving. In the
typical, textbook situation, you have an equation, say
3x + 2 = y
and you are asked to solve the equation for x. This

involves rearranging the symbols of the equation in the
familiar ways, e.g., moving the 2 to the right hand side
and dividing by the 3. These steps, originally termed
balancing and reduction are summarized in the
original meaning of the arabic word al-jabr3 used by 3
That is, .
Muhammad ibn Musa al-Khowarizmi (c. 780-850) in his
Compendious Book on Calculation by Completion and
Balancing. This is where our word algebra originates.
High school students are also taught a variety of ad
hoc techniques for solving in particular situations. For
example, the quadratic equation ax2 + bx + c = 0 can be
solved by application of the procedures of factoring, or
completing the square, or use of the quadratic formula:
p
b b2 4ac
x= .
2a
Parts of this formula can be traced back to at least the
year 628 in the writings of Brahmagupta, an Indian
mathematician, but the complete formula seems to date
from Simon Stevin in Europe in 1594 and was published
by Ren Descartes in 1637.
For some problems, students are taught named
operations that involve the inverse of functions. For
instance, to solve sin( x ) = y, one simply writes down
x = arcsin(y) without any detail on how to find arcsin

beyond use a calculator or, in the old days, use a table
from a book.
From Equations to Zeros of Functions

With all of this emphasis on procedures such as factoring
and moving symbols back and forth around an = sign,
students naturally ask, How do I solve equations in R?
The answer is surprisingly simple, but to understand
it, you need to have a different perspective on what it
means to solve and where the concept of equation
comes in.
The general form of the problem that is typically used
in numerical calculations on the computer is that the
equation to be solved is really a function to be inverted.
That is, for numerical computation, the problem should
be stated like this:
You have a function f ( x ). You happen to know the form of the
function f and the value of the output y for some unknown
input value x. Your task: find the input x that will produce
output y.
One way to solve such problems is to find an expression

for the inverse of f . The expression for the inverse can be
difficult or impossible to derive, but the notion of the
inverse function is simple. Theres even a standard
notation 4 to indicate the inverse of a function: f 1 . 4
Many students understand-
Fortunately, its rarely necessary to find an expression for ably but mistakenly take f 1
to mean 1/ f ( x ). The first is a
the inverse. Instead, the problem can be handled by
function inverse, the second
finding the zeros of f . is the value of the function
If you can plot out the function f ( x ) for a range of x, divided into 1.
you can easily find the zeros. Just find where the x where
the function crosses the y-axis. This works for any
function, even ones that are so complicated that there
arent algebraic procedures for finding a solution.
To illustrate, consider the function
plotFun(sin(x^2)*cos(sqrt(x^4+3)-x^2)-x+1 ~ x, 4
x.lim=range(-3,3)) 2
You can see easily enough that the function crosses the 0
y axis somewhere between x = 1 and x = 2. You can get -2
-2 -1 0 1 2
x
solving 43
more detail by zooming in around the approximate

solution:
0.5
plotFun(sin(x^2)*(cos(sqrt(x^4+3)-x^2))-x+1 ~ x, 0.0
x.lim=range(1,2))
-0.5
The crossing is at roughly x 1.6. You could, of -1.0
course, zoom in further to get a better approximation. Or, -1.5
you can let the software do this for you:

1.2 1.4 1.6 1.8
x
findZeros(sin(x^2)*(cos(sqrt(x^4+3)-x^2))-x+1 ~ x,
x.lim=range(1,2))
x
1 1.558
The syntax of findZeros() is very much like

plotFun(), but now the argument x.lim is used to state
where to look for a solution.
You need only have a rough idea of where the
solution is. For example:
x.lim=range(-1000,1000))
x
1 1.558
You can even say, I have no idea at all, by telling the

software to look anywhere between and .
x.lim=range(-Inf,Inf))
x
1 1.558
The findZeros() function will never look outside the

interval you specify. It will do a more precise job within
the interval if you can state the interval in a narrow way.
Setting up a Problem
As the name suggests, findZeros() finds the zeros of
functions. You can set up any solution problem in this
form. For example, suppose you want to solve ekt = 2bt

for b, knowing k, for example, k = 0.00035. You may, of
course, remember how to do this problem using
logarithms. But heres the set up for findZeros():
findZeros( exp(k*t) - 2^(b*t) ~ b, k=0.00035, t=1,

b.lim=range(-Inf,Inf) )
b
1 5e-04
Note that the I have no idea interval of to

was used. Youre usually better off if you have a finite
interval in mind, which you can set after you get a rough
idea of the solution:
findZeros( exp(k*t) - 2^(b*t) ~ b, k=0.00035, t=1,

b.lim=range(0,0.001) )
b
1 5e-04
Multiple Solutions
The findZeros() function will try to find multiple
solutions if they exist. For instance, the equation
sin( x ) = 0.35 has an infinite number of solutions. Here
are some of them:
findZeros( sin(x)-0.35 ~ x, x.lim=range(-10,10) )
x
1 -12.2088
2 -9.7824
3 -5.9256
4 -3.4992
5 0.3576
6 2.7840
7 6.6408
8 9.0672
solving 45
Exercises
Problem 1 (b) 4x2 2x + 20

Solve the equation A x = 3.33 or 1
sin(cos( x2 ) x ) x = 0.5 for x. B x = 3.33 or 1
0.0000 0.1328 0.2098 0.3654 0.4217 C x = 3.33 or 1
D x = 3.33 or 1
Problem 2 E No zeros
Find any zeros of the function
3e t/5 sin( 2p
2 t ) that are between t = 1 and (c) 2x3 4x2 3x 10
t = 10. Which one of these is a zero?
A There arent any zeros in that -1.0627 0 1.5432 1.8011 2.1223 3.0363 non
interval.
(d) 7x4 2x3 4x2 3x 10
B There arent any zeros at all! Which one of these is a zero?
C 2, 4, 6, 8 -1.0627 0 1.5432 1.8011 2.1223 3.0363 non
D 1, 3, 5, 7, 9
(e) 6x5 7x4 2x3 4x2 3x 10
E 1, 2, 3, 4, 5, 6, 7, 8, 9
Which one of these is a zero?
-1.0627 0 1.5432 1.8011 2.1223 3.0363 non
Problem 3
Use findZeros() to find the zeros of each Exercise 4 Construct a smoother for the
of these polynomials: height of pine trees as a function of age:
(a) 3x2 + 7x 10
pine = fetchData("Loblolly")
A x = 3.33 or 1
heightfun = spliner(height ~ age, data=pine
B x = 3.33 or 1
C x = 3.33 or 1 Use findZeros() to find the age at which
D x = 3.33 or 1 the height will be 35 feet.
E No zeros 11.3 11.9 12.2 12.7 13.1
4.2 Linear Algebra and Projection

Linear algebra operations are among the most important
in science and technology:
Project a single vector onto the space defined by a set A linear combination means to
of vectors. scale and add up vectors. The
matching problem consists of
Match a given vector with a linear combination of finding the correct scales so that
other vectors. the sum equals the given vector.
Projecting a vector onto a
In performing these operations, you will use two main space is analogous to seeing
functions, project() and mat(), along with the ordinary where the shadow of an object
will fall. It means to find the
multiplication * and addition + operations. There is also a point in the space that is closest
new sort of operation that provides a compact description to the vector.
for taking a linear combination: matrix multiplication,
written %*%.
To start, consider the sort of linear algebra problem
often presented in textbooks in the form of simultaneous
linear equations. For example:
x + 5y = 1
2x + 2y = 1 .
4x + 0y = 1
Many people will think of the above as a system of

three simultaneous equations. Another perspective is
valuable, however: treat the system as a single equation
involving vector quantities. To highlight the vectors,
re-write the equation like this:
0 1 0 1 0 1
1 5 1
x@ 2 A+y@ 2 A = @ 1 A.
4 0 1
Solving0 this vector

1 equation involves projecting the
1
vector ~b = @ 1 A onto the space defined by the two
01 1 0 1
1 5
vectors ~v1 = @ 2 A and ~v2 = @ 2 A. The solution, x
4 0
and y will be the number of multiples of their respective
vectors needed to reach the projected vectors.
solving 47
When setting this up with the R notation that you will

be using, you need to create each of the vectors ~b, ~v1 , and
~v2 . Heres how:
b = c(1,1,1)
v1 = c(1,2,4)
v2 = c(5,-2,0)
The projection is accomplished using the project()

function:
project(b ~ v1 + v2)
[,1]
v1 0.32895
v2 0.09211
Read this as project ~b onto the subspace defined by

~v1 and ~v1 .
The output is given in the form of the multiplier on ~v1
and ~v2 , that is, the values of x and y in the original
problem. This answer is the best in the sense that these
particular values for x and y are the ones that come the
closest to ~b, that is, the linear combination that give the
projection of ~b onto the subspace defined by ~v1 and ~v2 .
If you want to see what that projection is, just
multiply the coefficients by the vectors and add them up.
In other words, take the linear combination
0.32894737*v1 + 0.09210526*v2
[1] 0.7895 0.4737 1.3158
Notice that the projected value is not exactly

the same as ~b, even though it is the linear combination of
~v1 and ~v2 that reaches as close to ~b as possible, in other
words the projection of ~b onto the subspace spanned by
~v1 and ~v2 . The difference between ~b and its projection is
called the residual. The residual is also a vector and can
be calculated by subtracting the projection from ~b:
b - (0.32894737*v1 + 0.09210526*v2)
[1] 0.2105 0.5263 -0.3158
Matrices: Collections of Vectors

When there are lots of vectors involved in the linear
combination, its easier to be able to refer to all of them
by a single object name. The mat() function takes the
vectors and packages them together into a matrix. It
works just like project(), but doesnt involve the vector
thats being projected onto the subspace. Like this:
A = mat( ~ v1 + v2)
A
v1 v2
[1,] 1 5
[2,] 2 -2
[3,] 4 0
Notice that A doesnt have any new information; its

just the two vectors ~v1 and ~v2 placed side by side.
Heres the same projection, but using ~v1 and ~v2
packaged up together in A:
x = project( b ~ A)
x
[,1]
Av1 0.32895
Av2 0.09211
To get the linear combination of the vectors in A, you

matrix-multiply the matrix A times the solution x. The
odd-looking %*% operator invokes matrix multiplication.
A %*% x
[,1]
[1,] 0.7895
[2,] 0.4737
[3,] 1.3158
solving 49
The output from the matrix multiplication is, of

course, the same answer you got when you did the
vector-wise multiplication by hand.
The Intercept
Very often, your projections will involve a vector of all 1s.
This vector is so common that it has a name, the
"intercept." There is even a special notation for the
intercept in the mat() and project() functions: +1. For
instance:
A = mat( ~ v1 + v2 + 1)
A
(Intercept) v1 v2
[1,] 1 1 5
[2,] 1 2 -2
[3,] 1 4 0
Redundancy and the best solution

The instruction to project is implicit in traditional linear
algebraic notation. Rather than an imperative command
to do something, the traditional notation writes out the
relationship, e.g. A ~x = ~b, and leaves it to the human
reader to determine what quantities are unknown and
what operation is appropriate to find the unknown
quantities. For instance, if you know A and ~x, then
finding the unknown ~b in A ~x = ~b means performing
matrix multiplication A %*% x.
When, as in the example introducing this chapter, its
~x thats unknown, it often happens that that there is no
exact solution or that there is no unique solution. The
method used by project() looks for a best solution in
such cases. The appropriate meaning of best can
depend on the purposes for which the unknown is being
sought.
Part of the definition of best thats implicit in the
solution method used by project() is to make the
residual vector as small as possible: the so-called
least-squares solution. This is usually the standard and

appropriate thing to do.
Another aspect of best thats not so standard has to
do with handling cases where the vectors on the
right-hand side of ~ are linearly dependent, that is,
where there are more vectors than strictly needed to
define a unique subspace. As a simple example, consider
this projection problem:
b = c(3,5,-1)
v1 = c(1,2,3)
v2 = c(2,4,6)
The vectors v1 and v2 point in the same direction.

That is, even though there are two vectors v1 and v2, the
two vectors define a subspace that is only one
dimensional.
The project() method of solution handles this
redundancy by involving both vectors in the linear The statistical convention is
combination: implemented by the lm() op-
erator which works much like
project(b~v1+v2) project(). In addition to ex-
cluding redundant vectors, lm()
[,1] includes an intercept vector
unless explicitly forbidden to
v1 0.1429
do so by a -1 term, as in this
v2 0.2857
projection of ~b onto ~v1 and ~v2 :
Another possible choice, commonly used in statistics, lm(b~v1+v2-1)

is to consider the minimal set of vectors that define the
subspace and to exclude the other vectors from the
v1 v2
projection.
0.7143 NA
The exclusion of a redundant

Functions from Fitting vector is signalled by returning
NA as its coefficient in the linear
Often, model functions are created by fitting a linear combination.
combination of simple functions to data. The coefficients The outputs from project()
of such a fitted linear combination can be calculated using and from lm() are equivalent:
they imply exactly the same
project().
location for the projection of
Sometimes its convenient to package up the results of ~b, but the statistical method
a fit as a function rather than a vector of coefficients. To is more compact in that it
illustrate, consider the data contained in the involves fewer vectors in the
"cardata.csv" file giving measurements on various linear combination.
1978-79 model cars:
solving 51
cars = fetchData("cardata.csv")
Figure 12: Some of the data on

mpg pounds horsepower cylinders tons car size, power, and mileage
1 16.9 3968 155 8 2.0 from "cardata.csv"
2 15.5 3689 142 8 1.8
3 19.2 3281 125 8 1.6
4 18.5 3585 150 8 1.8
5 30.0 1961 68 4 1.0
6 27.5 2330 95 4 1.2
Suppose you want to construct a model function that

relates miles-per-gallon to the weight and horsepower of
the cars. You can, of course, use project():
project(mpg ~1+pounds+horsepower, data=cars)
[,1]
(Intercept) 46.932738
pounds -0.002902
horsepower -0.144931
The result is a vector of coefficients.

Compare the behavior of project() to that of
fitModel(). Whereas project() returns a vector of
coefficients, fitModel() gives back a function with named
parameters. The two are, of course, closely related. For
example:
fmpg = fitModel(mpg~A+B*pounds+C*horsepower, data=cars)

coef(fmpg)
A B C
46.932745 -0.002902 -0.144931
You may wonder, why use project, which works only

for linear parameters, when fitModel can handle both
linear parameters and nonlinear parameters.
Unfortunately, the mathematics of nonlinear fitting is
more difficult than that of linear algebra. This shows up
by the need to give initial guesses for the parameters in
fitModel, both linear and nonlinear. In contrast, the
linear algebra methods are simpler and do not require
any initial guess for the linear parameters. In addition,

models with only linear parameters are generally easier
to interpret and often provide a satisfactory
representation of the pattern seen in data. As such, its
very common to use linear models even if a nonlinear
model fits somewhat better. There is even an intermediate
sort of method, called generalized linear models (GLM)
that provides much of the benefit of linear models and
much of the flexibility of nonlinear models. Perhaps the
most commonly used GLM is called logistic regression.
Exercises
Exercise 1 YOUR TASK: Using the project() function:

Remember all those find the line that
goes through the points problems from
(a) Find the line that goes through the two
algebra class. They can be a bit simpler with
points (9, 1) and (3, 7).
the proper linear-algebra tools.
Example: Find the line that goes A y = x+2
through the points (2, 3) and (7, 8). B y = x + 10
One way to interpret this is that we are C y = x+0
looking for a relationship between x and y D y = x+0
such that y = mx + b. In vector terms, this E y =x 2
means that the x-coordinates of thetwo
2
points, 2 and 7, made into a vector (b) Find the line that goes through the origin
7
will (0, 0) and (2, 2).
be scaled by m, and an intercept vector
1
will be scaled by b. A y = x+2
1
B y = x + 10
x = c(2,7)
C y = x+0
y = c(3,-8)
D y = x+0
project(y ~ x + 1)
E y =x 2
[,1]
(Intercept) 7.4
x -2.2
(c) Find the line that goes through (1, 3) and
(7, 9).
Now you know m and b.
A y = x+2
B y = x + 10
C y = x+0
D y = x+0
E y =x 2
solving 53
Using project(), solve these sets of
(a) Find x, y, and z that solve the following:
simultaneous linear equations for x, y, and z:
0 1 0 1 0 1 0 1
1 5 1 1 1. Two equations in two unknowns:
x@ 2 A+y@ 2 A+z@ 2 A = @ 1 A.
4 0 3 1 x + 2y = 1
3x + 2y = 7
Whats the value of x?
-0.2353 0.1617 0.4264 1.3235 1.5739 A x = 3 and y = 1
B x = 1 and y = 3
(b) Find x, y, and z that solve the following: C x = 3 and y = 3
0 1 0 1 0 1 0 1 2. Three equations in three unknowns:
1 5 1 1
x@ 2 A+y@ 2 A+z@ 2 A = @ 4 A. x + 2y + 7z = 1
4 0 3 3 3x + 2y + 2z = 7
2x + 3y + z = 7
Whats the value of x?
-0.2353 0.1617 0.4264 1.3235 1.5739 A x = 3.1644, y = 0.8767, z =
0.8082
B x = 0.8767,y = 0.8082, z =
3.1644
C x = 0.8082, y = 3.1644, z =
0.8767
3. Four equations in four unknowns:
x + 2y + 7z + 8w = 1
3x + 2y + 2z + 2w = 7
2x + 3y + z + w = 7
x + 5y + 3z + w = 3
A x = 5.500, y = 7.356, z =
3.6918, w = 1.1096
B x = 1.1096, y = 3.6918, z =
7.356, w = 5.500
C x = 5.500, y = 7.356, z =
1.1096, w = 3.6918
D x = 1.1096, y = 7.356, z =
5.500, w = 3.6918
4. Three equations in four unknowns:
x + 2y + 7z + 8w = 1
3x + 2y + 2z + 2w = 7
2x + 3y + z + w = 7
True or False There is an exact solution.
(Hint: Whats the residual?)
5. Derivatives & Differentiation
The operator for taking derivatives is D(). Like all

ooperators, D() takes inputs and produces an output. D()
is quite simple: it takes just one input.
Input: an expression using the ~ notation. Examples:
x^2~x sin(x^2) y*cos(x)~y
On the left of the ~ is a mathematical expression written

in correct R notation that will evaluate to a number when
numerical values are available for all of the quantities
referenced. On the right of the ~ is the variable with
respect to which the derivative is to be taken. By no
means need this be called x or y; any valid variable name
is allowed.
The output produced by D() is a function. The
function will list as arguments all of the variables
contained in the input expression. You can then evaluate
the output function for particular numerical values of the
arguments in order to find the value of the derivative
function.
For example:
g = D(x^2 ~ x) g(1) g(3.5)
[1] 2 [1] 7
Formulas and Numerical Difference. When the expression

is relatively simple and composed of basic mathematical
functions, D() will often return a function that contains a
mathematical formula. For instance, in the above example
derivatives 55
function (x)
2 * x
For other input expressions, D() will return a function

that is based on a numerical approximation to the
derivative you cant see the derivative, but it is there
inside the numerical approximation method:
h = D( sin( abs(x-3) ) ~ x )
h
function (x)
numerical.first.partial(.function, .wrt, .hstep, match.call())
<environment: 0x1035a86a8>
Symbolic Parameters. You can include symbolic

parameters in an expression being input to D(), for
example:
s2 = D( A*sin(2*pi*t/P) + C ~ t)
The s2( ) function thus created will work like any

other mathematical function, but you will need to specify
numerical values for the symbolic parameters when you
evaluate the function:
s2
function (t, A, P, C)
A * (cos(2 * pi * t/P) * (2 * pi/P))
s2( t=3, A=2, P=10, C=4 )

s2(t, A = 2, P = 10, C = 4)
[1] -0.3883
1.0
plotFun( s2(t,A=2,P=10,C=4) ~ t, t.lim=range(0,20)) 0.5
0.0
-0.5
Partial Derivatives. The derivatives computed by D() are -1.0
partial derivatives. That is, they are derivatives where the

5 10 15
variable on the right-hand side of ~ is changed and all t
other variables are held constant.
Second Derivatives. A second derivative is merely the

derivative of a derivative. second derivatives You can use
the D() operator twice to find a second derivative, like
this.
df = D( sin(x) ~ x )
ddf = D( df(x) ~ x )
To save typing, particularly when there is more than

one variable involved in the expression, you can put
multiple variables to the right of the ~ sign, as in this
second derivative with respect to x:
another.ddf = D( sin(x) ~ x & x )
This form for second and higher-order derivatives also

delivers more accurate computations.
derivatives 57
Exercises
Exercise 1 Exercise 4 What will be the value of this

Using D(), find the derivative of derivative?
3*x^2 - 2*x + 4 ~ x. D( fred^2 ~ ginger )
(a) What is the value of the derivative at A 0 everywhere
x = 0? -6 -4 -3 -2 0 2 3 4 6 B 1 everywhere
C A positive sloping line
(b) What does a graph of the derivative
D A negative sloping line
function look like?
A A negative sloping line
Exercise 5
B A positive sloping line
Use D() to find the 3rd derivative of
C An upward-facing parabola
cos(2*t). If you do this by using the ~t&t&t
D A downward-facing parabola
notation, you will be able to read off a
formula for the 3rd derivative. What is it?
Exercise 2
A sin(t)
Using D(), find the derivative of
B sin(2t)
5*exp(.2*x) ~ x.
C 4 sin(2t)
(a) What is the value of the derivative at D 8 sin(2t)
x = 0? E 16 sin(2t)
-5 -2 -1 0 1 2 5 . Whats the 4th derivative?
A cos(t)
(b) Plot out both the original exponential B cos(2t)
expression and its derivative. How are C 4 cos(2t)
they related to each other? D 8 cos(2t)
A They are the same function E 16 cos(2t)
B Same exponential shape, but
different initial values
C The derivative has a faster Exercise 6
exponential increase Compute and graph the 4th derivative of
D The derivative shows an expo- cos(2*t^2)~t from t = 0 to 5. What does the
nential decay graph look like?
A A constant
Exercise 3 B A cosine whose period de-
2
Use D() to find the derivative of e x creases as t gets bigger
with respect to x (that is, exp(-(x^2) ~ x). C A cosine whose amplitude
Graph the derivative from x = 2 to 2. What increases and whose period
does the graph look like? decreases as t gets bigger
A A bell-shaped mountain D A cosine whose amplitude
B Exponential growth decreases and whose period
C A positive wave followed by a increases as t gets bigger
negative wave For cos(2*t^2)~t the fourth derivate is a
D A negative wave followed by a complicated-looking expression made up of
positive wave simpler expressions. What functions appear
in the complicated expression? derivative with respect to y, and these two

A sin and cos functions mixed partials:
B sos, squaring, multiplication pxy = D(x*sin(y) ~ x&y)
and addition pyx = D(x*sin(y) ~ y&x)
C cos, sin, squaring, multiplica-
tion and addition Pick several ( x, y) pairs and evaluate
D log, cos, sin, squaring, multipli- each of the derivative functions at them. Use
cation and addition the results to answer the following:
The partials with respect to x and to y are
identical. True or False
Exercise 7
Consider the expression x*sin(y) The second partials with respect to x and
involving variables x and y. Use D() to to y are identical. True or False
compute several derivative functions: the The two mixed partials are identical. That
partial with respect to x, the partial with is, it doesnt matter whether you
respect to y, the second partial derivative differentiate first with respect to x and
with respect to x, the second partial then y, or vice versa. True or False
6. Integrals & Anti-Differentiation
Youve already seen the D operator, which takes the

derivative of functions. Theres also an antiD operator,
which, as the name suggests, "undoes" differentiation. For
instance:
f = makeFun( x^2 ~ x )
Function Derivative Anti-Deriv. of Derivative

f df=D(f(x)~x) DF=antiD(df(x)~x)
4 4 4
3 2 3
DF(x)
df(x)
f(x)
2 0 2
1 -2 1
0 -4 0
-1 0 1 -1 0 1 -1 0 1
x x x
g = makeFun( 1/(1+exp(-x))~x )
Function Derivative Anti-Deriv. of Derivative

g dg=D(g(x)~x) DG=antiD(dg(x)~x)
1.0 0.25
0.4
0.8 0.20
0.2
DG(x)
0.6 0.15
dg(x)
g(x)
0.0
0.4 0.10
-0.2
0.2 0.05
-0.4
0.0 0.00
-4 -2 0 2 4 -4 -2 0 2 4 -4 -2 0 2 4
x x x
You can apply anti-differentiation to any function, not

just functions produced by differentiation. Indeed, its
rarely the case that you will want to anti-differentiate a

function that you have just differentiated. One undoes the
other, so there is little point except to illustrate how
differentiation and anti-differentiation are related to one
Flow (cups/sec)
another. But it often happens that the information you 0.10
have is a function that gives the derivative of some other

function of interest. 0.05
For example, consider how to measure the water lost in

a leak from a pipe. Its easy to measure the flow from the 0 10 20 30 40 50 60
Clock time (min)

leak; put a cup to catch the leaking fluid, then use a
stopwatch to time how long it takes to fill up the cup. Figure 13: The measured flow
Imagine that you have done this occasionally and from the pipe.
collected the following data over an hour:
time (min) 0 5 11 15 22 28 39 45 52 61
filltime (sec) 7 9 12 14 16 17 18 18 19 17
The flow rate has units of cups per second; you can
calculate it as 1/FillTime. Evidently, the leak is slowing at
4
first, then levels out, then seems to be increasing toward
Volume (cups/s min)

the end of the measurements. 3
You know the flow, but what is the total amount of 2
water that has leaked. Among the possibilities: 1
Instead of a cup, put a bucket under the leak and catch 0
all the water. Then you can measure how much is in 10 20 30 40 50
the bucket. Time (min)
Use the anti-derivative of flow. Figure 14: The accumulating

volume of the leak.
Remember, flow is the time derivative of volume.
Therefore, volume is the anti-derivative of flow with
respect to time.
To construct the anti-derivative function, you first
must put the derivative in the form of a function. A
spline through the measurements seems suitable. (The
data have been stored as "startR/leak-example.csv".)
leakdata = fetchData("startR/leak-example.csv")
flow = spliner( 1/filltime ~ time, data=leakdata)
volume = antiD( flow(time)~time )
The anti-derivative function thus calculated, volume, is

first and foremost a function of time. The derivative of
anti-derivatives 61
the function is the flow function, which is why the slope

of volume is steep at first and then climbs steadily.
Now suppose you wanted to know how much water
has leaked during a particular interval of time, say from
10 minutes to 30 minutes. Easy. Put a empty bucket
under the leak at time 10 minutes, then remove it at time
30 minutes and measure how much water has
accumulated. Or, leave the bucket where it is, mark how
much water is in the bucket at 10 minutes, record again at
30 minutes, and subtract. This second method is the
strategy by which an anti-derivative function is used to
find the accumulation:
In interpreting this number,
volume(time=30)-volume(time=10) 1.335, its important, as al-
ways, to take into account the
[1] 1.335 units. The output of an anti-
derivative function always
has units that are the prod-
Keep in mind that using the anti-derivative in this way
uct of the units of the deriva-
to compute an integral, accomplishes the same thing you tive function (cups/second
could do by accumulating the flow in a bucket, or in this example) and the vari-
multiplying the average flow by the duration. able of integration (minutes
Traditionally the definite integral is written like this here). So the number is re-
ally 1.335 cupminute/second.
(taking f (t) to be the flow): This may seem odd for a vol-
Z 30 ume, but remember that 60
f (t)dt. seconds equals a minute, so a
10 minute/second is just 60 sec-
ond/second is just 60. So the
There is also a notation for an indefinite integral: volume is 1.335 60 = 80.1
Z cups.
f (t)dt.
The two notations are very similar because the operations

are closely related. But keep in mind that the two things
are different.
The operation
R implied by writing the indefinite
integral f (t)dt is anti-differentation: the construction of
a new function. In the computer notation, constructing
this function is done with antiD.
f = flow # Just renaming to match the traditional notation

F = antiD( f(t) ~ t)
(Traditionally, capital letters are used to denote the

anti-derivative function.)
The operation implied by writing the definite integral

R 30
10 f ( t ) dt is implemented in three stages: construction of
the anti-derivative function, evaluation of the
anti-derivative function at the limits of integration, and
then taking the difference between these values. Those
last two stages are represented in this computer
command.
F(30) - F(10)
[1] 1.335
Exercises
Exercise 1 Find the numerical value of The function being integrated can have
each of the following definite integrals. additional variables or parameters beyond
R 5 1.5
1. x dx the variable of integration. To evaluate the
2
0.58 6.32 20.10 27.29 53.60 107.9 1486.8 definite integral, you need to specify values
R 10 for those additional variables.
2. 0
sin( x2 )dx For example, a very important function
0.58 6.32 20.10 27.29 53.60 107.9 1486.8 in statistics and physics is the Gaussian,
R 4 2x which has a bell-shaped graph:
3. 1
e dx
0.58 6.32 20.10 27.29 53.60 107.9 1486.8 gauss=makeFun((1/sqrt(2*pi*s^2))*exp(-(x-m)^2/
R 2 2| x | m=2,s=1.5)
4. 2
e dx plotFun(gauss(x)~x, x.lim=range(-10,10))
0.58 6.32 20.10 27.29 53.60 107.9 1486.8
0.25
Exercise 2 0.20
Theres a very simple relationship

gauss(x)
0.15
Rb Ra
between a f ( x )dx and b f ( x )dx 0.10
0.05
integrating the same function f , but 0.00
reversing the values of from and to. Create -5 0 5
some functions, integrate them, and x
experiment with them to find the (You might want to cut-and-paste this
relationship. definition of f() into your R session.) As you
A They are the same value. can see, its a function of x, but also of the
B One is twice the value of the parameters mean and sigma.
other. Evaluate each of the following definite
C One is negative the other. integrals:
D One is the square of the other. R1
1. 0 gauss( x, m = 0, s = 1)dx
Exercise 3 0.13 0.34 0.48 0.50 0.75 1.00
anti-derivatives 63
R2
2. 0
gauss( x, m = 0, s = 1)dx mathematical is represented as -Inf
0.13 0.34 0.48 0.50 0.75 1.00 on the computer.)
R2 0.13 0.34 0.48 0.50 0.75 1.00
3. 0 gauss( x, m = 0, s = 2)dx
R
0.13 0.34 0.48 0.50 0.75 1.00 5. gauss( x, m = 3, s = 10)dx

R3
4.
gauss( x, m = 3, s = 10)dx. (Hint: The 0.13 0.34 0.48 0.50 0.75 1.00
7. Dynamics
One of the most useful applications for the calculus

strategy involves studying how systems change in time:
dynamics. Often, dynamics are described in terms of
"differential equations." These relate the instantaneous
state of a system x (t) to the instantaneous change of state
dx
dt . Solving or "integrating" a differential equation
amounts to finding x (t) given dx dt .
As with other forms of integration, the setting is that
you have information about the derivative of a function
and want to deduce from this the function itself. When
the information given is a function of the independent
variable, for instance dx dt = f ( t ), then anti-differentiation
can be used. But the usual setting for differential
equations is that the derivative is known in terms of the
state itself:
dx
= g ( x ).
dt
The integrateODE() function solves such differential
equations. In order to complete the solution, the starting,
initial value of the state must be given. 0
To illustrate the difference between solving differential
1
equations and finding anti-derivatives, consider these two
x
different systems: 2
Anti-Derivative: The derivative is known as a function 4
of t. 1 2 3
dx Time t
= f (t) = 1 t.
dt Figure 15: Solution to dx
=1 t
dt
To find the solution, use antiD() on f (t) :
F = antiD( 1-t ~ t)
dynamics 65
Differential Equation: The derivative is known as a

function of x.
dx
= g( x ) = 1 x.
dt
To find the solution, use integrateODE() and specify how
the change in the state x depends on x. 1.0
0.8
soln = integrateODE( dx ~ 1-x, x=0,
State x
0.6
tdur=list(from=0,to=4)) 0.4
0.2
Note that its necessary to give an initial condition for 0.0
the state; the argument x=0 sets the initial condition to 1 2 3
zero. You also have to specify the time interval over Time t
which the solution is to be found.
Figure 16: Solution to dx
dt =
Although the two settings involve similar-looking
1 x from initial condition
formulas, they lead to completely different results. x=0
Example: Logistic growth. The dynamical model for

logistic growth involves a state x that grows at small
values, but decays at large values:
dx
= rx (1 x/K ).
dt
For x < K the population grows while for x > K the
population decays. The state x = K is a stable
equilibrium." Its an equilbrium because, when x = K, the
change of state is nil: dx/dt = 0.
The integrateODE() function takes the differential
equation as an input, together with the initial value of the
state. Numerical values for all parameters must be
specified, as they would in any case to draw a graph of
the solution. In addition, must specify the range of time
for which you want the function x (t). For example, heres
the solution for time running from 0 to 20.
soln <- integrateODE( dx ~ r*x*(1-x/K),

x=1, K=10, r=.5,
tdur=list(from=0,to=20))
The object that is created by integrateODE() is a

function of time. Or, rather, it is a set of solutions, one for
each of the state variables. In the logistic equation, there
is only one state variable x. Finding the value of x at time
t means evaluating the function at some value of t. Here

are the values at t = 0, 1, . . . , 5:
soln$x(0:5)
[1] 1.000 1.548 2.320 3.324 4.509 5.751

10
Often, you will plot out the solution against time: 8
soln$x(t)
6
plotFun( soln$x(t)~t, t.lim=range(0,20))
4
Differential equation systems with more than one

state variable can be handled as well. To illustrate, here is 5 10 15
t
the SIR model of the spread of epidemics, in which the
state is the number of susceptibles S and the number of Figure 17: A solution to
infectives I in the population. Susceptibles become logistic-growth dynamics.
infective by meeting an infective, infectives recover and
leave the system. There is one equation for the change in
S and a corresponding equation for the change in I. The
initial I = 1, corresponding to the start of the epidemic.
epi = integrateODE( dS~-a*S*I, dI~a*S*I - b*I,

a=0.0026,b=.5,S=762,I=1,tdur=20)
This system of differential equations is solved to

produce two functions, S(t) and I (t).
600
plotFun( epi$S(t)~t, t.lim=range(0,20))
epi$S(t)
plotFun( epi$I(t)~t, add=TRUE, col="red") 400
200
In the solution, you can see the epidemic grow to a
peak near t = 5. At this point, the number of susceptibles 0
has fallen so sharply that the number of infectives starts 5 10 15
t
to fall as well. In the end, almost every susceptible has
been infected. Figure 18: A solution S(t) and
I (t) to the SIR dynamics.
The phase plane provides a powerful format to
visualize how the state of a two-variable system changes
in time. Each point in the phase plane represents a state
of the system. The dynamical rule the differential
equations describe a flow field. The evolving state of
the system is a trajectory that starts at an initial point and
follows the flow.
dynamics 67
The graphics for phase-plane display are best

implemented using a graphical user interface. One such
interface, is provided by the mPP() program (which works
with RStudio). First, you must load the mPP software.
fetchData("mPP.R")
fetchData("DiffEQ.R")
Differential equations are specified using a function.

For instance, the SIR equations can be implemented like
this:
SIR = function(suscept,infective){
a=0.0026; b=0.5
dsuscept = -a*suscept*infective;
dinfective = a*suscept*infective - b*infective;
return( c(dsuscept, dinfective) );
}
To run the graphical user interface, call mPP():
mPP(DE=SIR, xlim=c(0,1000), ylim=c(-0,1000))
Figure 19: The mPP display with

the SIR model.
Example: Another Dive from the Board Consider a

diver as she jumps off the high board and plunges into 5
the water. You construct a dynamical model with state

4
Height (m)
3
variables v and x: velocity and height. 2
1
dive = integrateODE( dv~-9.8, dx~v, 0
v=1,x=5,tdur=1.2) 1
plotFun( dive$x(t)~t, t.lim=range(0,1.2)) 0.2 0.4 0.6 0.8 1.0
t
Whats nice about the differential equation format is
Figure 20: Diving without
that its easy to add features like the buoyancy of water buoyancy or drag.
and drag of the water. Well do that here by changing the
acceleration (the dv term) so that when x < 0 the
acceleration is slightly positive with a drag term
proportional to v2 in the direction opposed to the motion. 4
Height (m)
2
diveFloat = integrateODE(
dv~ifelse( h>0, -9.8, 1-sign(v)*v^2), dh~v, 0
v=1,h=5,tdur=10) 2
plotFun( diveFloat$h(t)~t, t.lim=range(0,10))
2 4 6 8
t
According to the model, the diver resurfaces after
slightly more than 5 seconds, and then bobs in the water. Figure 21: Buoyancy returns
the diver to the surface, drag
damps out the oscillations.
Exercises
Exercise 1 Use integrateODE() to find a of transmission. Suppose you could half

solution to dx
dt = 0.3x with initial condition the rate of transmission, by, say, closing
x = 2. Find the solution value at time t = 10. schools. What happens to the peak
number of infectives?
Exercise 2 Consider the SIR dynamics
epi=integrateODE(dS~ -a*S*I, dI~a*S*I-b*I,
a=0.0026,b=.5,S=762,I=1,tdur=20) Exercise 3 A simple model for the principal
P of a home mortgage at 5% annual interest
1. Instead of starting with 1 infective (I=1), is dP
dt = 0.05P 10000 where the annual
there are 10 infectives. Find the peak payment is on the mortgage is $10000 and t
number of infectives as the epidemic is in years. Starting from an initial value of
plays out. P = 100000, find the time at which the
2. The parameter a corresponds to the rate mortgage value goes to zero.
Activities
Activity 1: Where is the landscape the steepest?

Youre going to create a terrain that you can examine
visually. Your goal is to calculate where the terrain is
steepest.
1. Create a function representing the height H of a terrain
and plot it out in contour form. By changing the seed,
you can create your own terrain.
H = rfun(~x&y, seed=677) 4 -4
plotFun(H(x=x,y=y)~x&y,
2
x.lim=range(-5,5),y.lim=range(-5,5))
-2
-2 0
0 2
0
y
2 6
-2
2. The partial derivatives reflect the slope of the terrain in 8
the cardinal directions along the axes. You can -4 4
calculate each of the partials H/x and H/y. For -4 -2 0 2 4
example x
dHdx = D(H(x=x,y=y)~x)
3. Plot out H/x as a function of x and y.

4
plotFun(dHdx(x=x,y=y)~x&y, 4
2
x.lim=range(-5,5),y.lim=range(-5,5)) -4 0 2
0
y
0 -2
2
-2
4. Steepness can be either uphill or downhill,
-2
0
depending on the direction you are heading. To -4
measure steepness, as opposed to slope, you might -4 -2 0 2 4
sensibly take the absolute value of the slope or the x

square of the slope. Plot out the absolute value
|H/x | as a function of x and y. (Hint: the abs()
function.) Plot out also the square: (H/x )2
5. The above plots reflect the steepness only in the x

direction. To get overall steepness, add up the
steepness in the x and in the y directions. You can do
this
p as either |H/x | + |H/y| or as
(H/x )2 + (H/y)2 .
6. Looking at your plot of overall steepness, find the
steepest point in the original terrain, that is, the x and
y at which the terrain is steepest.
For discussion:
Why the square-root in the expression for steepness?
(Hint: Think about the units.)
Which do you think is more appropriate for measuring
steepness: the sum of absolute values of slopes or the
square-root of the sum of square slopes?
Is it hard to find the x and y at which the function are
steepest? How does the smoothness of the terrain
compare to the smoothness of the steepness function?
Question for a computer programming class: Given a

terrain function and initial location x0 and y0 as inputs,
find the path followed by a drop of water as it rolls down
the landscape.
Activity 2: How far does the car get?

Suppose you want to calculate the capacity of a road
intersection regulated by a traffic light. One part of this
calculation involves the distance that a car can travel,
starting from a standstill, in a given amount of time.
Much of the information you have about the car is in
the form of velocity. Its common sense that the velocity is
0 mph at the start and you can reasonably presume that
the car accelerates to the speed limit of the road, say 30
mph. By doing a little observation, you can get an
estimate for how many seconds it takes a typical car to
reach the speed limit. For the present, lets say that its 6
seconds.
Construct a function (of time t) that models the
velocity of the car and is consistent with the
activities 71
information specified above. Implement this as a

computer function and graph it to verify that the
velocity profile meets the specification. Here are
several possible forms (with parameters that might be
inappropriate).
v1 = makeFun(pmin(30, 30*t/6) ~ t)
v2 = makeFun(30*pnorm(t,mean=3,sd=1)~t)
v3 = makeFun(ifelse(t>6,30,30*t/6) ~ t)
Compute the integral of the velocity function to get

position. Plot out the position as a function of time.
(Question: What should the from argument be? The to
argument?)
Does the result how far the car goes depend
strongly on the details of the velocity function? Is there
a strong reason to prefer one form of the function to
the others? Make a quantitative argument either way.
Activity 3: Allocating resources for health care

The data file "jmm2012data1.csv" contains experts
evaluation of the number of Quality-Adjusted Life Years
(QALYs) that will result from investment in two different
health care options, A and B.
You can read in the file like this:
dat = fetchData("jmm2012data1.csv")
The data are discrete. You can construct a continuous,

smooth function of QALYs versus expenditure for each of
A and B like this:
fA = spliner(A~expend,data=dat,monotonic=TRUE)
fB = spliner(B~expend,data=dat,monotonic=TRUE)
You can now use fA and fB like any other function, for
example, plotting
150
plotFun(fA(xA)~xA, xA.lim=range(0,50))
fA(xA)
100
The problem. You have a total budget of 50 units to 50
spend between A and B.

0
10 20 30 40
xA
(a) Whats the best allocation of the funding between A

and B to maximize the QALY output?
(b) If your budget were increased slightly, what would be
the resulting change in optimal QALY output? Put
this in the form of a fraction: budget units per
increased QALY.
(c) Suppose that you were mandated to spend at least 40
units of the budget on A. How would that additional
constraint change the resulting optimal QALY output?
Activity 4
Here is an AP Calculus exam problem published by the
College Board. p
If the function f is defined by f ( x ) = x3 + 2, and g
is an antiderivative of f such that g(3) = 5, then g(1) =?
(A) 3.268
(B) 1.585
(C) 1.732
(D) 6.585
(E) 11.585
You can use the antiD() operator to compute the
antiderivative.
g = antiD(sqrt(x^2+2) ~ x)
The antiderivative is usually treated as an expression

that involves some unknown additive constant C. One
way to answer the question is to find the value of C
consistent with g(3) = 5, then apply this in calculating
g (1).
Instructor Notes
This document is for instructors: people who are already

very familiar with the mathematical concepts of calculus
but who are new to R or perhaps generally to the use of
computer software in teaching calculus.
Realize at the outset that computer notation you will
be using is somewhat different from traditional
mathematical notation:
Operators are often written differently. Examples: x^2
or sqrt(x) or A*sin(2*pi*t/5). Students will make
many mistakes at first, for example writing 2 t instead
of 2*t, or leaving out parentheses to contain
arguments (e.g., writing sin x rather than sin(x)).
You give names to results so that you can use them in
new places.
Variable names are often multi-letter, latin, and use
words rather than diacritical marks. Example: instead
of x or p0 , write perhaps xestim or popInitial.
An explicit functional style is used where the operator
name is followed by a pair of parentheses. Within the
parentheses, the various arguments are separated by
commas, and often identified byname with a notation
like x=3.
This last point marks an important difference between the
computer notation and traditional notation, since in
traditional position rather than a label is used to identify
different inputs. For instance, in traditional notation you
Rb
can write a f ( x )dx. In this notation, the operator is
integration, a and b are arguments that mark the lower
and upper bounds of integration, and depending on your
philosophy, either f ( x ) is the quantity being integrated
(with dx denoting that x is the variable of integration), or

f ( x )dx is the input to the integration operation. Similarly,
df
in dx , the operation is differentiation, f is the function
being differentiated, and the variable with respect to
which the derivative is taken is identified by the dx as x.
In many calculus operations there is both a
mathematical function and one or more variables with
respect to which an operation is being performed. In
English we often will say with respect to, or as a
function of or versus. To mimic this in R, we have
adopted a built-in R-language capability called a formula.
The central component of a formula is the ~ symbol:
tilde.
To illustrate, here is the construction of a derivative:
D( sin(x^2) ~ x )
function (x)
cos(x^2) * (2 * x)
The computer has performed the derivative,

producing a function of x.
Typically, you would perform this operation in order
to use the resulting function for some purpose, so you
should give the output a name:
fp = D( sin(x^2) ~ x )
To evaluate the function at a particular value of x, use

the parentheses to pass along to fp the desired value for
x:
fp( x=2 )
[1] -2.615
Or, perhaps you want to plot out the function for

some range of x values:
5
plotFun( fp(x) ~ x, x.lim=range(0,4) )

fp(x)
You might reasonably conclude that ~ turns an

-5
expression into a function. Thats not quite right. The ~
provides a place to identify one or more symbolic
1 2 3
x
instructor notes 75
variables with respect to which the operation

differentiation, plotting, etc. is to take place. The ~
symbol is a way to organize symbolic information, not an
instruction to perform a particular calculus operation.
You must always specify explicitly the operation you
want to perform. If, for example, you want to create a
function, you can do so with the makeFun operator:
f = makeFun( sin(x^2) ~ x )
Notice that the operation has been named f, not f ( x )

or some such thing. When you want to use the function,
you must give it an argument. If the argument is
numerical, just go ahead and put the numerical value in
parentheses:
f(x=3)
[1] 0.4121
If the argument is a symbol, and youre planning to

use that symbol as the variable in some calculus
operation, then you will want to use the ~ marker as well,
for example:
plotFun( f(x) ~ x, x.lim=range(0,4) )
or
F = antiD( f(x) ~ x )
Multiple Variables
In modeling-based calculus, its appropriate to introduce
functions of multiple variables very early. The notation is
designed to make this straightforward. For example:
plotFun( exp(-t/10)*sin(2*pi*x/5) ~ x&t, 8

0.5 -0.5
x.lim=range(0,5), t.lim=range(0,10) ) 6
0.0
t
4
When you create a function of multiple variables,
theres a question of how to identify which variable is 2
which. One traditional approach is to use a notation that

1 2 3 4
x
conveys both the name of the function and the

arguments, e.g., f ( x, t) = e t/10 sin 2p
5 x.
In the computer notation, the name is just f . The
arguments are conveyed in the definition.
f = makeFun( exp(-t/10)*sin(2*pi*x/5) ~ x&t )
Perhaps it seems obvious from the definition that the

first argument will be x and the second t. That happens
to be true here, but its not something to rely on. Its
better to identify the arguments explicitly, e.g.
f( x=1, t=2 )
[1] 0.7787
f( t=2, x=1 )
[1] 0.7787
When you are using a function in a symbolic context,

this can lead to a notation that looks, at first glance, to be
0.6
redundant:
f(t = t, x = 2)
0.5
plotFun( f(t=t,x=2)~t, t.lim=range(0,10))
0.4
or, for example, to integrate with respect to x but leave 0.3

things as a function of symbolic t:
0.2
F = antiD( f(t=t,x=x)~x ) 2 4 6 8
Symbolic Parameters
You can use symbols as parameters. The makeFun(), D(),
antiD(), and plotFun() operators will recognize them
and keep them in symbolic form. For instance:
D( A*exp(-k*t) ~ t )
function (t, A, k)
-(A * (exp(-k * t) * k))
However, when you want to evaluate such a function,

including plotting it, you need to assign a numerical
value to the parameter.
instructor notes 77
gp = D( A*exp(-k*t) ~ t )
plotFun( gp(t,A=2,k=1/10) ~ t, t.lim=range(0,20))
gp(t, A = 2, k = 1/10)
-0.05
-0.10
The Basic Calculus Operations -0.15
Differentiation and Partial Differentiation -0.20
5 10 15
t
D( a*x^2 + b*x + c ~ x )
function (x, a, b, c)
a * (2 * x) + b
D( a*x^2 + b*x + c ~ x&x )
function (x, a, b, c)
a * 2
D( a*x + b*y + c*x*y ~ x&y )
function (x, y, a, b, c)
c
Anti-Differentiation/Indefinite Integration
F = antiD(exp(x)*x^2 ~ x)
Evaluating an Anti-Derivative/Definite Integration

The anti-derivative function produced by antiD() has the
same arguments as the input expression.
To compute an integral, evaluate the anti-derivative at
the top and bottom of the interval of integration and
subtract:
F(x=2) - F(x=1)
[1] 12.06
When integrating to , use the values -Inf or Inf.

The output of antiD() is arranged so that the

"constant of integration" appears as an argument with a
default value of zero. only the upper bound appears in
the argument list. When a numerical approach is taken,
the lower bound will be set to a default value. You can set
the lower bound explicitly, if you wish, for instance:
F = antiD(exp(x)*x^2 ~ x)
F(x=2) # default const. of integration is 0
[1] 12.78
F(x=2, C=10) # setting const. of integration to 10
[1] 22.78
For some expressions, antiD() will compute the

anti-derivative in symbolic form; the function returned
will have familiar R notation. For all other input
expressions, the output of antiD() is still a function, but
it will not have an algebraic form. Instead, the function
will compute an integral numerically.
The integral has an numerical interface that makes it
operate as if it were an indefinite integral, but in fact the
underlying calculation is always definite integration. In
particularly, there is a default lower bound of zero. You
need be concerned with the lower bound only if the
default of zero is inappropriate. This can occur if there is
a singularity in the expression being anti-differentiated at
or near zero. You can specify a different lower bound by
using the argument lower.bound when calling antiD().
You can use symbolic parameters in defining an
anti-derivative, but they must be given numerical values
at the time of evaluating it.
Fbell = antiD( dnorm(x,m=m,s=s) ~ x ) # symbolic s and m

Fbell(x=12,s=2,m=10) - Fbell(x=8,s=2,m=10)
[1] 0.6827
instructor notes 79
2.3 makeFun versus function

If you have used R extensively before, you may wonder
why makeFun() is being used rather than the more
conventional function(). In all the documentation for R,
after all, the standard way to define a function is
something like this:
f = function(x,y) { A + B*x^2 + C*y^2}
This is a very useful format for computer

programmers. It provides an explicit list of the
arguments, sets the order of the arguments, and allows
the functions body to refer to parameters that are defined
in an enclosing environment, e.g. the global environment.
The corresponding use of makeFun() goes like this:
g = makeFun(A + B*x^2 + C*y^2 ~ x&y)
The arguments are being identified as the variables

following the ~. But when makeFun() packages up a
function, it parses the body of the function for all variable
names and includes them as formal, named arguments to The exception is the name pi,
the function: which is treated as a constant
with value p.
g
function (x, y, A, B, C)
A + B * x^2 + C * y^2
This has some advantages for teaching mathematics

and statistics.
Students can successfully use the system without
having to learn or even be aware of the scoping rules
of the language.
The notation reinforces that variables and parameters
are not intrinsically different.
It avoids giving special status to an operator that
happens to be called function and allows students to
see as natural other operators that produce functions
as output, e.g., fitModel(), D(), antiD(), rfun(), and
so on.
Index
abs, 69 findZeros, 4345 loblolly pine, 28

algebra, 41 fitModel, 34, 37, 38, 51, 79 logistic growth, 65
anti-derivative, 64 fitting, 34 logistic regression, 52
anti-differentation, 77 flow vs volume, 60
anti-differentiation, 59 function, 79 makeFun, 16, 18, 79
antiD, 64, 72, 7779 function, 79 makeFun, 79
axis labels, 23 mat, 46, 48, 49
axis ranges, 23 generalized linear model, matrix, 46, 48
52 max, 32
car prices, model, 36 graphing functions, 16 min, 32
coef, 35 monotonic spline, 32
coefficients, model, 34 half-life, 37 mosaic, 8
completing the square, 41 mosaic package, 16
constant of integration, 77 Income-Housing data, 20, mPP, 67
contour plots, 25 36 Muhammad ibn Musa
cubic spline, 28 indefinite integral, 61 al-Khowarizmi, 41
installing packages, 16
D, 5458, 79 integrateODE, 64, 65, 68 names, 21
data files, 20 integration, 59 Newton, Isaac, 4
definite integral, 61 constant of, 77 non-parametric function, 29
derivative, 54 numerical, 78 nonlinear parameter, 34
Descartes, 41 interaction, 36 NSF, 8
differential equations, 64 intercept term, 49 numerical calculus, 5
differentiation, 54, 77 interpolating function, 28 numerical integration, 78
diving, model, 68 inverse functions, 42
dynamics, 64 parameter
lHopital, Marquis de, 4 nonlinear, 34
epidemics, 66 least-squares solution, 50 parameter vs variable, 79
equations, simultaneous, 46 linear algebra, 46, 51 parameters, 16
exponential function, 37 linear combination, 46 partial derivatives, 55
linear dependence, 50 partial differentiation, 77
f, 62 linear function, 34 phase plane, 66
factoring, 41 linear interpolant, 29 plotFun, 16, 18, 19, 22, 23,
fetchData, 20 lm, 50 25, 43
index 81
plotPoints, 23 simultaneous linear tides, hawaii, 40

plots equations, 46
contour, 25 SIR model, 66
units
surface, 25 smoother
of integral, 61
plotting functions, 16 function, 28
Used Honda data, 35
project, 4653 span, 30, 31
Utilities data, 34
projection, 46, 48 solving equations, 41
spline
quadratic equation, 41 cubic, 28 variable vs parameter, 79
function, 28 variables, in datasets, 21
integrating, 60 vector, 46
redundancy of vectors, 49 monotonic, 32
residuals, 38 spreadsheets, 20
rfun, 79 Wagon, Stan, 40
subspace, 48
RStudio, 11, 67 water, cooling, 38
sum of square residuals, 38
surface plots, 25
scatterplots, 20 symbolic calculus, 5 zero-crossing, 43
scoping rules, 79 symbolic parameters, 55, 76 zero-finding, 41

Start R in Calculus

Uploaded by

Document Informationclick to expand document informationLearning calculus with the R programming language

Document Informationclick to expand document information

Copyright:

Available Formats

Start R in Calculus

Uploaded by

Document Information

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Start R in Calculus

Uploaded by

Copyright:

Available Formats

DANIEL KAPLAN

Copyright (c) 2013 by Daniel Kaplan.

Cover Photo: . [Photo credit: Maya Hanna.]

To Netta, Liat, and Tamar.

Its no coincidence that calculus and calculate are

forme cette regie.

produit des autres.

adx. Celle de a -t- x x b

These rules are still useful in important ways, but there

on aura xcalculus, 6c comme ces deux

understand and carry out. As such, it provides a means

R is a language for communicating instructions to a

There are several aspects of technical comput-

not the thought process that you went through in com-

If you have ever travelled to a country where you dont

or Chinese or Spanish, you will make mistakes and you

My thanks go to my colleagues at Macalester, par-

Functions & Graphing 16

R is a language that you can use to direct a computer to

Figure 2: The RStudio window

The RStudio window has the familiar menu bar and is

commands. The console maintains a transcript of your

When you press return, R interprets your command

Simple arithmetic in R is done with a familiar

5-2 3^2 9*6+1

[1] 3 [1] 9 [1] 55

5/2 9*7 9*(6+1)

[1] 2.5 [1] 63 [1] 63

For functions with letter names, for instance cos and

sqrt(9) acos(0) log(10)

[1] 3 [1] 1.571 [1] 2.303

cos(1) sin(pi) exp(-3)

[1] 0.5403 [1] 1.225e-16 [1] 0.04979

Figure 4: Use parentheses to

When you want to see the value of a named value, just

Of course, R can do much more than such

Important Note about Typing: In the RStudio

Exercise 1 Use R to calculate a numerical Exercise 3 Here is a set of R expressions

2.1 Graphing Mathematical Functions

2. Real-world quantities are not typically named x and y,

the bounds of the horizontal axis over which you want

The plotFun() operator puts this all together, taking 25

Often, its natural to write such relationships with the 2 4 6 8

Sometimes, you want to give a function a name so that

Once the function is named, you can evaluate it by

To plot g(), you need to provide an argument and use

Exercise 1 This command Exercise 4

plotFun(). Hand in the command that you

(a) 4x 7 in the window x from 0 to 10.

2.2 Making Scatterplots from Data

Income IncomePercentile CrimeProblem AbandonedBuildings

All of the variables in the dataset will be shown

[1] "Income" "IncomePercentile" "CrimeProblem"

When you want to access one of the variables, you

[1] 3914 10817 21097 34548 51941 72079

[1] 39.6 32.4 26.7 23.9 21.4 19.9

Even though the output from names() shows the

Spelling and capitalization are important. If you make

5/2 97 9(6+1)

plotFun( sin(2pit/10)exp(-.2x) ~ t & x,

g = makeFun( sin(2pit/10)exp(-.2x) ~ t & x) 6