\documentclass[11pt]{article}
\def\Z{\mathbb{Z}}
\def\Q{\mathbb{Q}}
\usepackage{latexsym,amsmath,amssymb,epsfig}
\setlength{\textwidth}{6.5in}
\setlength{\oddsidemargin}{0in}
\setlength{\textheight}{8.0in}
\begin{document}
\centerline{MONT 105N -- Analyzing Environmental Data}
\centerline{Chapter 10 Project}
\centerline{March 18, 2019}
\vskip 10pt
\noindent
\subsection*{Introduction}
The goals of this chapter's project are:
\begin{itemize}
\item To introduce some basic features of the {\tt R} statistics 
package.  This is a much more extensive and powerful piece of software
than the statistical routines in Google Sheets or Excel. In fact {\tt R} is 
probably the go-to ``industry standard'' software for academic statisticians at the moment. 
\item To see how {\tt R} computes some descriptive statistics 
and graphics for understanding the ``shape'' or ``distribution'' 
of a data set. 
\item To gain an intuitive understanding of the \emph{Central Limit Theorem}, 
a basic result in statistics that underlies the wide applicability of normally-distributed
random variables.  
\end{itemize}

\subsection*{Background on {\tt R} -- Getting Started}
The instructions here are geared toward the Windows version
of {\tt R}.   {\tt R} is a free,
open source package and you can obtain a copy to run
on your own computer, and there are Windows, 
Mac-OS, and Linux versions available.  These have the same 
functionality but slightly different interfaces.   

On PC's when you double-click the {\tt R} shortcut from the
standard desktop, you will launch the R-GUI (graphical user interface) window.
Inside this is a workspace window with a white background.
It \emph{only} allows you to 
\begin{itemize}
\item enter commands, and 
\item view numerical output.
\end{itemize}
It does not let you create integrated documents, and
graphics are displayed in separate windows within the 
R-GUI, as you will see.

You can \emph{print} the workspace window or any graphics
window at any time.  Just highlight that window and
press the toolbar printer button.  

To save your work, with comments and answers to questions
(e.g. for the project report), I recommend 
\begin{itemize}
\item opening another editor window (either Word or NotePad
are fine for this on Windows systems), and
\item using it to create a text file that will be
your permanent record of your work.
\item When you have a result you want to save from
the {\tt R} workspace, copy both the input command and
the output and paste into the text file.
\item Any time you want to add a comment to explain
a calculation or answer a question, enter that 
directly into the text file. 
\item Save that text file print it when you want a hardcopy
(from Word or NotePad).  
\end{itemize}

\subsection*{A First {\tt R} Session}
\vskip 10pt
In this first session, you will basically work through some 
examples illustrating features of {\tt R} and some of the 
descriptive statistics we have talked about so far in 
the text.  Keep a record of what you do for future reference.

In its most basic form, {\tt R} provides a command-driven
``statistical calculator'' in which you type in data and
commands to generate numerical and graphical output that
are displayed.  When you have the {\tt R} window open, note the
$>$ input prompt after the opening message.  This is 
where you should start. 

(1)  The following data set records
systolic blood pressure readings from 15 adults (not a very
healthy group, on the whole!)
First we need to get the systolic blood pressure readings
into R so that we can work with them.  The basic way is to 
enter (after the $>$ prompt):
\vskip 10pt
\begin{center}
{\tt bpdata <- c(172,140,123,130,155,148,138,129,137,161,123,152,133,128,142)}
\end{center}
\noindent
and press ENTER.  Note:  the $<-$ is formed by typing the $<$ and $-$
characters right next to each other.  It is the {\it assignment operator}
in R.  This takes the data on the right and assigns it to the name on
the left.   The $c$, followed by the list
of numbers, separated by commas, is {\tt R}'s way of constructing 
{\it an ordered list or vector}.  You should just get another
input prompt if everything went OK with this (error messages if not).
If you did get an error, try again, and try to follow the above exactly.
You can verify that the numbers are entered correctly by typing in the 
name {\tt bpdata} as a command; the numbers will be displayed in a row.

(a)  We can compute the mean $\overline{y}$ and standard deviation 
$s$ of the data in {\tt R} (in several different ways!):  Type in 
\vskip 10pt
\begin{center}
{\tt ybar <- sum(bpdata)/length(bpdata)}\\
{\tt mean(bpdata)}\\
{\tt s <- sqrt(sum((bpdata - ybar)\^{}2)/(length(bpdata) - 1))}\\
{\tt sd(bpdata)}
\end{center}
\noindent
and press ENTER after each line to see the results.  Note the way the sum function
works in a very nice way here---in the computation of {\tt ybar}, we only need
to put in the name of the data list; the sum of all the elements is computed by 
default.  Then in the computation of {\tt s}, we can subtract \emph{each entry}
of the {\tt bpdata} vector from {\tt ybar} as shown!  Of course we will usually
want now to use the built-in functions {\tt mean} and {\tt sd}, but
the first of each pair of computations shows how we can also use more primitive functions
in {\tt R} when we need them.  The command {\tt summary(bpdata)} automates the computation 
of the five number summary if you want to see that.
\index{five number summary}

(b)  Chebyshev's theorem mentioned in Chapter 9 is a result about the distribution of data 
sets.  It says in a particular case
that \emph{there will always be at least $75\%$ of the numbers in a data set 
within 2 standard deviations of the mean}.  Let's use {\tt R} to check that this is true here.  In its most basic
use, {\tt R} can function as a simple numerical calculator.
\vskip 10pt
\begin{center}
{\tt upperb <- ybar + 2*s}\\
{\tt lowerb <- ybar - 2*s}
\end{center}
\noindent
How many of the data values are in the range $[\overline{y} - 2 s, \overline{y} + 2s]$,
though?  We could just count by hand, of course, with a data set this small.
But {\tt R} also contains a very rich collection of facilities for selecting subsets of a 
data set and massaging data in various ways.  Here is a basic way to select the portion 
of the data set in the range we want:
\vskip 10pt
\begin{center}
{\tt middle <- subset(bpdata, bpdata <= upperb \& bpdata  >= lowerb)}\\
{\tt length(middle)/length(bpdata)}
\end{center}
(The first line picks out the subset of the data satisfying both of the inequalities; the second
computes the fraction of the whole data set contained in the ``middle.'')
Did Chebyshev's theorem ``work'' here?
\index{Chebyshev's theorem}


(2)  Now let's look at a different data set.  This gives the number of parasites found in 
each of a sample of 100 individual foxes.
The information is given there in table format:
\begin{center}
\begin{tabular}{l|ccccccccc}
Parasites & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8\\ \hline
Foxes & 69 & 17 & 6 & 3 & 1 & 2 & 1 & 0 & 1\\
\end{tabular}
\end{center}
So we want to create a data set containing 69 0 entries, 17 1 entries, etc.
Clearly this is going to be tedious and error prone if we do it the most basic
way(!)  Fortunately, {\tt R} lets us construct lists in more flexible ways
besides just listing the entries.  We can indicate how many times to repeat
a given value if we want, like this:
\vskip 10pt
\begin{center}
{\tt foxpar <- rep(c(0,1,2,3,4,5,6,7,8),c(69,17,6,3,1,2,1,0,1))}\\
\end{center}
\noindent
The {\tt rep} stands for ``repetition'' here.
Note that the first list inside the {\tt rep} gives the values of the number of parasites---the first row 
in the table format---and the second list gives number of times each value appears in the data.  
The first list could also be abbreviated as
{\tt 0:8} ({\tt R}'s notation for a list of successive integer values).  

\index{frequency histogram!relative}
\index{histogram}
(a)  To construct a relative frequency histogram
for the number of parasites per fox in {\tt R}, let's begin by entering
\vskip 10pt
\begin{center}
{\tt hist(foxpar)}\\
\end{center}
\noindent
The histogram will be generated in a separate graphics window.  
This is a ``default'' version of the command.  The {\tt hist}
command takes a number of options that control how the bins
are selected and how the histogram is drawn as well.  You 
will often want to modify the default results.   For
instance, looking at the output, how were the bins chosen here?
Is that optimal for showing the different numbers of foxes
with each number of parasites?  Suppose we want bins centered
on the integer values $0,1,2,3,4,5,6,7,8$ as we did in other
examples.  We can put in an option to make
{\tt R} do it that way too.  The {\tt breaks} option takes
a list of values and makes those the boundaries of the bins 
for the histogram:
\vskip 10pt
\begin{center}
{\tt hist(foxpar,breaks=c(-0.5,0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5))}\\
\end{center}
\noindent
How is this different and why is it better than the default
way that {\tt R} drew the histogram?  

(b) To practice using the commands we talked about
before, calculate $\overline{y}$ and $s$ for this data set.
How many of the numbers were within two standard deviations
of the mean in this case?  

(c)  The list of breaks for the histogram can be
anything in {\tt R}.  What happens if you change the command above to 
\vskip 10pt
\begin{center}
{\tt hist(foxpar,breaks=c(-0.5,0.5,1.5,2.5,3.5,4.5,5.5,8.5))}\\
\end{center}
\noindent
(deleting the 6.5 and 7.5 break points, so we're lumping
together the counts for the values 6,7,8 parasites)? 
What is different about the new histogram that {\tt R} drew in this case?  
Look at the graph carefully: do the heights of the boxes still represent 
\emph{frequencies}?  (Hint: what
they do represent are \emph{densities} that are computed
by taking the frequency count and dividing by the 
width of the box in each case.  What is true if we 
add the areas of the boxes drawn this way?  We will see
other reasons why this is a preferable way to draw histograms
with unequal bin sizes later. )
If we \emph{really want frequencies}, we can make {\tt R} do it that way:
\vskip 10pt
\begin{center}
{\tt hist(foxpar,breaks=c(-0.5,0.5,1.5,2.5,3.5,4.5,5.5,8.5),freq=TRUE)}\\
\end{center}
\noindent
However note that now {\tt R} generates a warning message saying
the AREAS in the plot are ``wrong.''  Think about this.
Why is the last box in this version of the histogram
\emph{misleading} if drawn this way?  

\subsection*{Background on Probability Distributions in {\tt R}}
The R package contains built-in functions for computing
probability functions for the discrete random 
variables we have discussed, plus functions for the probability
density functions and cumulative distributions of the 
standard types of continuous random variables.  Here's how it works. 
Each type of random 
variable is covered by a \emph{family} of 4 functions distinguished by a 
\emph{prefix letter} {\tt d, p, q}, or {\tt r}:
\vskip 10pt
\begin{itemize}
\item {\tt d} -- the p.m.f.  in the discrete case, or the p.d.f. in the continuous case,
\item {\tt p} -- the c.d.f.
\item {\tt q} -- the quantile function (it computes quartiles, deciles, percentiles; we'll see examples
of how it is used later)
\item {\tt r} -- random number generator (for sampling random values
of a random variable with that distribution)
\end{itemize}
Following the prefix letter comes the rest of the function 
name and the inputs needed to compute the corresponding values.  
\vskip 10pt
\noindent
\subsection*{Distributions of Sample Means}
\index{Central Limit Theorem|(}
Work through the following example before and then follow the 
pattern here in the other parts of question A below. 
You do not need to hand in any work for this preliminary
discussion, just for the parts of the question itself.
\begin{enumerate}
\item Begin by generating a list of $100$ random numbers from a particular
distribution.\footnote{This is a type of continuous random
variable that we did not discuss in the main part of the text.  It is one of the exponential
distributions from an example in the optional section \S 10.8.  You will see roughly what
its p.d.f. looks like below even if you did not look at the optional section.}  
Use the command:
\vskip 10pt
\begin{center}
{\tt x <- rexp(100,1/3)}
\end{center}
\noindent
to assign the result to the name {\tt x}.  We will think of 
this as a \emph{sample} of values from that distribution.
\item Plot the density histogram for this data set {\tt x} using the 
command
\vskip 10pt
\begin{center}
{\tt hist(x,freq=FALSE)}
\end{center}
This should look roughly like the part of the graph of the 
p.d.f.  for positive $y$, reflecting the fact that the random numbers
are generated according to probabilities given by that density function.
\item Our main goal is to understand what happens if we generate \emph{lots}
of such samples, find their sample means, then consider \emph{the distributions
of those sample means}.  
\item Here is one way to do this in {\tt R} using
a simple \emph{for loop}.  First we create space to store the means of 
1000 different random samples:
\vskip 10pt
\begin{center}
{\tt means <- array(1:1000)}
\end{center}
Next, we generate 1000 different random samples, find the 
means of each of them, and store them in the array created above:
\begin{verbatim}
                for (i in 1:1000) 
                   { 
                   x <- rexp(100,1/3)
                   means[i] <- mean(x)
                   } 
\end{verbatim}
(Note: This could all be entered as one input line provided you leave
spaces in the appropriate places.  However, I think it is more readable 
(and it is easier to identify typos if you happen to make one) 
if you press ENTER at the end of each of the lines as
you are typing this.  If you do it that way, you should note that {\tt R} generates a
new input prompt {\tt +} each time, indicating that you 
are still in the body of the {\tt for} loop.  The final {\tt \}}
will close off the loop and take you back to the $>$ prompt.)
\item What does the distribution of the sample means look like?
We can see that by generating a frequency histogram for the data in the {\tt means}
array:
\vskip 10pt
\begin{center}
{\tt hist(means,breaks=20,freq=FALSE)}
\end{center}
This should look completely different from the histogram showing
the distribution of the single sample above.  In fact, what
does this histogram remind you of?
\item If you said, ``a normal p.d.f.,'' good.  If you 
said something different, go back and look again.
\item Now, which normal density is it?  Well,  by Theorem 10.36
from the text, if 
we have 100 sampled values $Y_1, \ldots, Y_{100}$ from any distribution
where $E(Y_i) = \mu$ for all $i$, and $V(Y_i) = \sigma^2$ for all $i$, and
$$\overline{Y} = \frac{1}{100}(Y_1 + \cdots + Y_{100})$$
is the sample mean, then we will have
$$E(\overline{Y}) = \mu$$
and 
$$V(\overline{Y}) = \frac{\sigma^2}{100},$$
so the SD of the sample means will be $\sigma/10$.   
For the exponential random variables we are looking at, we will take it 
on faith that $\mu = 3$ and 
$\sigma = 3$.  Hence the expected value of the sample means is also 
$3$ and the SD of the sample means is $3/10 = .3$.
\item  Let's overlay that normal density and see how well things match:
\vskip 10pt
\begin{center}
{\tt curve(dnorm(x,3,.3),add=TRUE)}
\end{center}
(The {\tt dnorm} is the normal p.d.f. The 3 is the $\mu$ and the .3 is the $\sigma$.  You
will use other normal densities in the questions below.  They are specified in the same way.)
This should match quite well (not perfectly, though, since there will always be
variation in the sampled values as we know from the project for Chapter 9).
\end{enumerate}

\subsection*{Assignment} 
Repeat the parts of the worked example below for the means of 1000 random
samples of size 100 from each of the following distributions.  Note:  In
each case to find the appropriate normal density curve, you will need to 
recall facts about expected values and variances of the $Y_i$ individually,
then use Theorem 10.36.
\vskip 10pt
\begin{enumerate}
\item[(A)] A uniform distribution with $min = 4$ and $max = 10$.  (These
are the endpoints of the interval where the uniform density is different
from zero.)   Use {\tt runif(100,4,10)} to generate a random sample
of size 100 from this uniform distribution.
See Theorem 10.33 in the text for the appropriate $\mu$ and $\sigma^2$ in this case.
Then apply Theorem 10.36 as in the example above.  
\item[(B)]  A Poisson distribution with $\lambda = 2$.  From the discussion in the text, we know
$\mu = \lambda$ and $\sigma^2 = \lambda$ for Poisson random variables.
The appropriate {\tt R} function to generate the samples of size 100 is 
{\tt rpois(100,2)}.  
\end{enumerate}
What is happening in each case?  Explain.  Is the distribution of the sample
means approximately normal every time?  
\end{document}
