---
title: "Normal models"
output:
  html_document: default
  pdf_document: default
date: "9/1/2017"
---

```{r setup, include=FALSE}
library(knitr)
knitr::opts_chunk$set(echo = TRUE)
```

```{r top_pic, echo=F}
dx = 0.01
x <- seq(-3.5,3.5, by=dx)
y <- dnorm(x)
plot(x,y, type='l', ylim=c(0,0.6), yaxs='i', xlim=c(-3,3), xlab="", ylab="", yaxt='n', lwd=3)
x <- seq(-1,1, by=dx)
y <- dnorm(x)
polygon(c(x,rev(x)), c(rep(0, length(y)), rev(y)), col='lightgray')
```

Most statistical data is *modeled* by some *distribution*. By far the most common such distribution that arises in practice is the [normal distribution](https://en.wikipedia.org/wiki/Normal_distribution).

## Normal curves

### The standard normal

The standard normal distribution is a specific bell shaped curve; its graph is shown at the top of the page above. It has a maximum as it crosses the $y$-axis; its *mean* is zero. It decreases as we move away from the mean and tapers off to the $x$-axis. It's scaled so that the total area under its graph is one. The area under the curve is distributed in a specific way; it's *standard deviation* is one.

### More normals

The standard normal is one of a specific family of bell shaped curves. We can change the mean by shifting the graph to the left or to the right. We can change the standard deviation by compressing or dilating along the $x$-axis while simultaneously doing the opposite along the $y$-axes. Several such curves are shown in the figure below.

```{r normal_curve_comparisons, echo=F}
x <- seq(-5,5, length=400)
y <- dnorm(x, -2,0.4)
plot(x,y, 'l', lwd=2, col='blue', xlab="", ylab="")
y <- dnorm(x, 0,1)
lines(x,y, 'l', lwd=2, col='red')
y <- dnorm(x, 2,2)
lines(x,y, 'l', lwd=2, col='green')
text(-3.02,0.7, 'm=-2', col='blue')
text(-3,0.62, 's=0.4', col='blue')
text(0,0.55, 'm=0', col='red')
text(0,0.48, 's=1', col='red')
text(3.2,0.3, 'm=2', col='green')
text(3.2,0.22, 's=2', col='green')
```

Note that the mean and standard deviation are denoted $m$ and $s$ in the figure above. Often, we denote these with
the Greek letters $\mu$ and $\sigma$.

### A formula

It's worth mentioning that there is a specific formula that generates the normal curves, namely
$$f_{\mu,\sigma}(x) = \frac{1}{\sqrt{2\pi}\sigma} e^{-(x-\mu)^2/(2\sigma^2)}.$$
You can explore the graphs of these functions [on Desmos](https://www.desmos.com/calculator/cqkqhoe91k).

### Modeling histograms with normal curves

Often, if we scale a histogram so that the total area of its rectangles is one, then the normal curve with the same mean and standard deviation will match the histogram fairly closely. An example is shown below.

```{r, echo=F}
set.seed(2)
d = rnorm(1000,10,5)
hist(d, main='A (scaled) histogram and its normal model', xlab = '', ylab = '', col='gray', freq = F)
m = mean(d)
s = sd(d)

dx = 0.01
x <- seq(-10,30, by=dx)
y <- dnorm(x,m,s)
lines(x,y, lwd=3)
```

We've got a [dynamic illustration on our webpage](https://marksmath.org/visualization/normal_model/).

## Standard deviation as a ruler

The mean and standard deviation can be used to standardize data. A principle tool to do so is the $z$-score.

### The heptathlon

I've got a CSV on my website that contains the results of the 2008 Olympic Heptathlon. Here are the first few rows:

```{r}
df = read.csv('https://marksmath.org/data/heptathlon.csv')
kable(head(df,3))
```

Note that Nataliia Dobrynska won the long jump by about half a meter over the average jump; Hyleas Fountain won the 200m by a second and a half faster than the mean. Which of these two results deserves more points?

One way to compare these results involving different units is to use the standard deviation as a ruler. To do so, we standardize. Given a data point $X$ chosen from data that is normally distributed with mean $\mu$ and standard deviation $\sigma$ we compute
$$Z = \frac{x-\mu}{\sigma}.$$
With this choice of $Z$, note that 
$$X = \mu + Z\sigma.$$
Thus, $Z$ measures exactly how many standard deviations $X$ is from the mean. This computation is sometimes called the $Z$-score.

Here's how we apply this to the 200m vs long jump question:

```{r}
# Z-score for the 200m
two_hundread_x = min(df$two_hundred)
(two_hundread_x-mean(df$two_hundred))/sd(df$two_hundred)
```


```{r}
# Z-score for the long jump
long_jump_x = max(df$long_jump)
(long_jump_x-mean(df$long_jump))/sd(df$long_jump)
```

It looks like the long jump is the more impressive result.

### SAT Example

The SAT is designed to have mean score of 500 with a standard deviation of 100. Let's suppose you score a 650.

 a. What is your $z$-score?
 b. What is your percentile?

*Solution*: For part (a), the $z$-score is simply
$$Z = \frac{650-500}{100} = 1.5.$$
We can answer part (b) by computing an area under the standard normal curve as shown below:

```{r, echo = F}
dx = 0.01
x = seq(-4,4, by=dx)
y = dnorm(x)
plot(x,y, type='l', ylim=c(0,0.6), yaxs='i', xlim=c(-3,3), xlab="", ylab="", yaxt='n', lwd=3)
x <- seq(-4,1.5, by=dx)
y <- dnorm(x)
polygon(c(x,rev(x)), c(rep(0, length(y)), rev(y)), col='lightgray')
text(-0.5,0.15,"A = 0.933")
```

This type of area is typically computed using either software or a table. We will learn how to do this in R soon but, for now, let's take a look at [this table](https://marksmath.org/classes/Fall2017Stat185/handouts/tables.pdf).

## The 68-95-99.7 rule

There is a statistical rule of thumb called [the 68-95-99.7 rule](https://en.wikipedia.org/wiki/68%E2%80%9395%E2%80%9399.7_rule) that states that

  * 68% of the population lies within 1 standard deviation of the mean,
  * 95% of the population lies within 2 standard deviations of the mean, and
  * 99.7% of the population lies within 3 standard deviations of the mean.

For this to work, the property being measured should be normally distributed. If so, then the rule of thumb follows from the fact that this property holds for the standard normal:

```{r percentiles, echo=F, message=F}
library(openintro)
data(COL)
X <- seq(-4, 4, 0.01)
Y <- dnorm(X)
plot(X, Y,
     type = 'n',
     axes = FALSE,
     xlim = c(-3.2, 3.2),
     ylim = c(0, 0.4),
     xlab = "", ylab = "")
abline(h = 0, col = COL[6])
at <- -3:3
labels <- expression(mu - 3 * sigma,
                     mu - 2 * sigma,
                     mu - sigma,
                     mu,
                     mu + sigma,
                     mu + 2 * sigma,
                     mu + 3 * sigma)
axis(1, at, labels)
for (i in 3:1) {
  these <- (i - 1 <= X & X <= i)
  polygon(c(i - 1, X[these], i),
          c(0, Y[these], 0),
          col = COL[i],
          border = COL[i])
  these <- (-i <= X & X <= -i + 1)
  polygon(c(-i, X[these], -i + 1),
          c(0, Y[these], 0),
          col = COL[i],
          border = COL[i])
}

# _____ Label 99.7 _____ #
arrows(-3, 0.03,
       3, 0.03,
       code = 3,
       col = '#444444',
       length = 0.15)
text(0, 0.02, '99.7%', pos = 3)

# _____ Label 95 _____ #
arrows(-2, 0.13,
       2, 0.13,
       code = 3,
       col = '#444444',
       length = 0.15)
text(0, 0.12, '95%', pos = 3)

# _____ Label 68 _____ #
arrows(-1, 0.23,
       1, 0.23,
       code = 3,
       col = '#444444',
       length = 0.15)
text(0, 0.22, '68%', pos = 3)

lines(X, Y, col = '#888888')
abline(h = 0, col = '#AAAAAA')
```


### Example

Let's suppose that the mean life expectancy of a cat is 14 years with a standard deviation of 2.5 years. Assuming that the cats' life spans are normally distributed, is it reasonable to expect a cat to live to 22 years old?

*Solution*: Well, the $z$-score for a 22 year old kitty cat would be
$$Z = \frac{22-14}{2.5} = 3.2.$$
As we know, only $0.3\%$ of cats live beyond a $z$-score of 3, so a 22 year old cat would be quite rare indeed.

## Assessing normality

In the previous problem, we *assumed* that the lifespan of cats was normally distributed. Is this a valid assumption? I have no idea, *but* given a data set there are ways to check.

### A histogram

Human heights are classically known to be normal. Let's take a look at the heights of the over 9500 men chosen from our [CDC data set](https://marksmath.org/data/cdc.csv).

```{r, echo=F}
cdc = read.csv("https://marksmath.org/data/cdc.csv")
men = subset(cdc, gender=='m')
hist(men$height, col='gray', xlab='', ylab='', main='Histogram of male heights', breaks = 30)
```

Looks kinda normal. 

### The normal probability plot

There's a more sensitive tool, though: the *normal probability plot* - more generally called a *quantile-quantile plot*.

Here’s the basic idea: Suppose we have some sample data consisting of $n$
numerical values. We’ll make a plot with $n$ points - one point for each of the $n$ values.
The vertical component for a given data point is just the value of the data point itself. The horizontal component for that given data point is determined by the $z$-score of that percentile.

For example, $z=1$ is greater than about 84% of all values produced by the normal distribution. You can look that up in a normal table or use the following R command:

```{r}
pnorm(1)
```

Furthermore, the $84^{\text{th}}$ quantile for the men's CDC height data is 73 inches, as the following computation reveals.

```{r}
quantile(men$height,0.84)
```
Thus, the point $(1,73)$ should be on the normal probability plot for this data. Here's the normal probability plot of that height data, together with the point $(1,73)$ shown in yellow.

```{r, echo=F}
qqnorm(men$height)
qqline(men$height)
points(1,73, pch=21,  bg="yellow")
```

By contrast, income is classically not normally distributed. Here's a normal probability plot for a random sample of incomes taken from the US census.

```{r, echo=F}
cen = read.delim('https://www.marksmath.org/data/census.txt')
qqnorm(cen$totalFamilyIncome)
qqline(cen$totalFamilyIncome)
```
