\documentclass[11pt]{article}
\def\Z{\mathbb{Z}}
\def\Q{\mathbb{Q}}
\def\R{\mathbb{R}}
\def\C{\mathbb{C}}
\usepackage{latexsym,amsmath,amssymb}
\usepackage[all]{xy}
\setlength{\textwidth}{6.5in}
\setlength{\oddsidemargin}{0in}
\setlength{\textheight}{8.0in}
\begin{document}
\centerline{Mathematics 375 -- Probability Theory}
\centerline{Computer Lab Day 2 -- Probability Distributions in R}
\centerline{October 19, 2011}
\vskip 10pt
\noindent
{\it Goals}
\vskip 10pt
The goals of today's lab are:
\begin{itemize}
\item to gain some more practice with basic features of the R statistics 
package, and 
\item to reinforce intuition about cumulative distribution
functions (CDFs), and probability density functions (PDFs) for
continuous random variables.
\end{itemize}
\vskip 10pt
\noindent
{\it Basic R Mechanics}
\vskip 10pt
Refer to the assignment sheet for the first lab day
if you need to refresh your memory about running 
R in Haberlin 136, about how to save and print your
results, etc.  The Lab 1 assignment is posted on 
the course homepage, so you can open it and look
at the instructions there as needed.  
\vskip 10pt
\noindent
{\it Background on Probability Distributions in R}
\vskip 10pt
The R package contains built-in functions for computing probabilities
and cumulative density functions for many of the discrete random 
variables we have discussed, plus functions for the probability
densities and cumulative distributions of many of the 
standard types of continuous random variables that we will study 
this semester and next.  Here's how it works.  (Also look at table
on 332 of Dalgaard \emph{Introductory Statistics with R}.) 
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 probability function in the discrete case, or the 
density function continuous case,
\item {\tt p} -- the cumulative distribution
\item {\tt q} -- the quantile function (essentially the \emph{inverse
function} of the cumulative distribution, but of course the cumulative
distribution can fail to be injective, so this needs to be taken with
 a grain of salt)
\item {\tt r} -- random number generator 
\end{itemize}
Following the prefix letter comes the rest of the function 
name and the inputs needed to compute the corresponding values.
For instance the family of functions for Poisson random variables 
(discrete) looks like this:
\begin{itemize}
\item {\tt dpois(y,lambda)} -- the probability function $P(Y = y)$ for 
$Y \sim Poisson(\lambda)$; 
the first input is the $y$ (which must be an integer), the second is
the parameter $\lambda$.
\item {\tt ppois(y,lambda)} -- the cumulative distribution $P(Y \le y)$
\item {\tt qpois(q,lambda)} -- computes the $y$ such that $P(Y \le y) = q$
\item {\tt rpois(N,lambda)} -- generates a list of $N$ random numbers following
a $Poisson(\lambda)$ distribution
\end{itemize}
The distributions that we will use today have these names (after the
prefix letter explained above):
\begin{itemize}
\item {\tt binom} -- binomial (need to supply a $y,q$, or $N$ as above, and 
the $n,p$ values)
\item {\tt norm} -- normal (need to supply a $y,q$, or $N$ as above, and values for {\tt mean} and {\tt sd} parameters)
\item {\tt exp} -- exponential (need to supply a $y,q$, or $N$ as above, and a value for the {\tt rate} parameter)
\end{itemize}
\vskip 10pt
\noindent
{\it Today's Lab}
\vskip 10pt
\noindent 
{\it Visualizing Probability (or Density) and Distribution Functions}
\vskip 10pt 
Work through the following examples before attacking the question A. 
below.  You do not need to hand in any work for this preliminary
discussion, just for the question.
\begin{enumerate}
\item The \emph{probability histogram} for a discrete 
random variable that we saw in Chapter 3 can be viewed as 
a sort of graph of the probability function.  This can be 
can be drawn in {\tt R} as follows.  Suppose for instance that 
$Y \sim Binomial(50,.3)$.  The $y$ for which $P(Y = y) \ne 0$ 
are $0 \le y \le 50$, so we begin by creating that sequence of
$y$-values, then plot the {\tt dbinom} function like this:
\vskip 10pt
\begin{center}
{\tt y <- seq(0,50,1)}\\
{\tt plot(y,dbinom(y,50,.3),type="h")}
\end{center}
\vskip 10pt
\noindent
This should all be relatively self-explanatory after you execute
the commands and see the graph.  (If you don't understand what the $y$
is here, just enter the command {\tt y} to print out what is 
in that variable.)  You want to think of the vertical
``pins'' as representing the boxes in the probability histogram.
The {\tt type = "h"} indicates this style of plotting.  What happens
if you change that to {\tt type = "l"} (that's an ell, not the 
digit 1), or leave out the type option all together?
\item  Now, to see the cumulative distribution function, enter 
\vskip 10pt
\begin{center}
{\tt plot(y,pbinom(y,50,.3),type="h")}
\end{center}
Note how this relates to the general properties we discussed last time.
The same types of commands will work for any of the functions described above. 
For instance, these commands
\vskip 10pt
\begin{center}
{\tt y <- seq(-4,4,0.1)}\\
{\tt plot(y,dnorm(y,0,1),type="l")}
\end{center}
will plot the graph of the probability density of a 
random variable with a normal distribution with mean 0 and sd 
(standard deviation) 1.  Of course,
we have not discussed these in class yet, but for now you 
can just take this graph as a sort of \emph{definition} for what that
type of continuous random variable looks like.  
\item In the normal
family of random variables, the mean can be any real number and
the sd can be any positive number.  You may want to experiment with 
changing those and the {\tt y} in the commands immediately 
above so that you get a good understanding of how the density
function changes.  
\item As you do this, you should notice that
it is necessary to change the set of values supplied for
plotting to get reasonable graphs.  If you want to have {\tt R}
choose the values automatically, you can use another command
called {\tt curve}.  For instance
\vskip 10pt
\begin{center}
{\tt curve(dnorm(x,0,1))}
\end{center}
\vskip 10pt
\noindent
will do essentially the same thing as the sequence of
two commands for the {\tt dnorm(y,0,1)} function above.
{\it Note:}  In the {\tt curve} command the ``variable'' in
the {\tt d}-function (or {\tt p}-function) always
has to be called {\tt x}, though.  Another nice feature
of the {\tt curve} command is that you can add an option
{\tt add = T} inside the parentheses if you want to overlay
the curve on the active plot.  See p. 42-43 in the R book for
a detailed example.  
\end{enumerate}
\vskip 10pt
\noindent
Question A.  
\vskip 10pt
The probability histogram for $Y \sim Binomial(50,.3)$
is actually very close to the graph of the density function
for a normal random variable with some mean and sd.  Thinking
about what we know about binomial random variables, determine
the appropriate mean and sd values, and generate a composite
plot showing the probability histogram overlaid with the
graph of the normal density.  In the output you submit, 
include the commands that generated the output, and an explanation
of how you determined which mean and sd values to use.  
\vskip 10pt
\noindent
{\it Sampling from Given Distributions}
\vskip 10pt
The {\tt r-} function in the family for each type of 
distribution as described above generates lists of random
numbers following the corresponding distribution. This
means, for instance, that for the continuous types, 
the probability that each number comes out in
a given interval is determined by the density function:
$$P(a \le Y \le b) = \int_a^b f(y)\ dy.$$
In this problem, we will generate two collections of 
1000 random numbers, one
from a normal distribution with mean 1 and sd 1, and the
other from an exponential distribution with rate 1.  We have
not discussed these yet formally either, but the exact formula
for the exponential density is one we used in an example in
class on Friday:
$$f(y) = \begin{cases}  e^{-y} & y > 0\\
                          0 & y\le 0.
         \end{cases}$$ 
\vskip 10pt
\noindent
Question B.  
\begin{enumerate}
\item[1.]  You can get an idea of what this all means by plotting the probability densities as in question A.  Generate these plots and write a short description of the differences between these two distributions.  
\item[2.]  Enter the commands
\vskip 10pt
\begin{center}
{\tt nr <- rnorm(1000,1,1)}\\
{\tt er <- rexp(1000,1)}
\end{center}
\vskip 10pt
\noindent
These generate two collections of 1000 random numbers, one
from an normal distribution with mean 1 and sd 1, and the
other from an exponential distribution with rate 1. 
Suppose you did not know which list had been 
produced which way.  Thinking about what we know about
normal distributions from the ``empirical rule,'' determine
a way to test which is which by considering the intervals
$\overline{y} \pm k s$ for $k = 1,2,3, ... $.  Carry out
your calculations for both lists {\tt er} and {\tt nr}
and discuss your results.  \emph{Note:}  If you reexecute the 
{\tt rnorm} and {\tt rexp} commands, you will be generating
different lists of random numbers each time.  So the 
values of the (sample) mean and (sample) sd will be changing
each time.  However, they should be quite similar and 
consistent enough that your conclusion does not change.  
\end{enumerate}
\vskip 10pt
\noindent
{\it Assignment}
\vskip 10pt
Lab reports containing input, output, and answers to the questions
posed above are due in class no later than Monday, October 24.  
If you do not finish during the hour today, you can return
to HA 136 any time it is not in use by another class.
 \end{document}
