---
title: "Week 8 Lecture 1"
output: html_document
---

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

### Probability and coin tosses, again

Let's review coin tosses. Here, we'll generate $100$ coin tosses where the weight of the coin is 0.7 (i.e., we expect the coin to come up Heads (1) 7 times out of 10).

```{r}
tosses <- rbinom(n = 100, size = 1, prob = 0.7)
tosses
```

We can comfirm that the coin indeed came up heads about 70% of the time:
```{r}
mean(tosses)
```

What is the probability mass function here?

| Tails | Heads |
|-------|-------|
|  0.3  |  0.7  |


### The binomial distribution

Let's now try something different: we'll toss the coin 10 times at a time, and record how many times it came up Heads. Then we'll repeat *that* number 100 times.

The result will be 100 records of the number of heads out of 10 tosses:

```{r}
rbinom(n = 100, size = 10, prob = 0.7)
```

Let's now make a histogram out of this. We'll use 50000 tosses so that the histogram looks "smoother"

```{r}
set.seed(0)
exp10 <- rbinom(n = 500000, size = 10, prob = 0.7)
exp10.df <- data.frame(n.heads = exp10)
ggplot(exp10.df) + geom_bar(mapping = aes(x=n.heads, y = ..prop..)) + scale_x_continuous(breaks = 1:10)
```

We can take this histogram as representing (approximately) the probability mass function of the Binomial distribution with probability $0.70$ and 10 trials.

In other words, here are the probabilities of obtaining $k$ heads for each $k$ between $0$ and $10$:

```{r}
binom.pmf.df <- data.frame(n.heads = 0:10,
                  binom.pmf = dbinom(x = 0:10, size = 10, prob = 0.7))
binom.pmf.df
```

We didn't need to look the numbers up on the histogram (and the histogram was an approximation anyway). R can compute the PMF using a built-in formula. But it's nice that those match up!

Let's look more closely at `dbinom` for a second. `dbinom` computes the probability that you'll get a particular number of heads with `size` trials:

```{r}
dbinom(x = 2, size = 10, prob = 0.7)
```

You can give `pbinom` several inputs by sending a vector as the `x`

```{r}
dbinom(x = c(2, 3), size = 10, prob = 0.7)
```

Here, we have 11 possible outcomes. All probabilities will always add up to 1. Think of it this way: the probability of 0 or 1 or 2 or ... or 10 heads is 1 (since one of those is sure to happen). To get that probability, we could add up the probabilities of the individual events. In other words:

$P(\text{nH = 0 or nH = 1 or ... or nH = 10}) = P(\text{nH = 0}) + P(\text{nH = 1}) + ... + P(\text{nH = 10}) = 1$.

Here's another application of this idea. What is the probability of observing a coin coming up heads once or more?



We have
$$P(\text{nHeads} \geq 1) + P(\text{nHeads} = 0) = 1$$

since one of those events will definitely happen and the events are mutually exclusive. So we have

$$P(\text{nHeads} \geq 1)  = 1 - P(\text{nHeads} = 0)$$

and

```{r}
1 - dbinom(x = 0, size = 10, prob = 0.7)
```

We can use the same idea to compute the probability that the coin will come up heads more than once:

 $\text{P(nH = 2)} + ... + \text{P(nH = 10)} = 1 - \text{P(nH = 0)} - \text{P(nH = 1)}$


```{r}
1 - dbinom(x = 0, size = 10, prob = 0.7) - dbinom(x = 1, size = 10, prob = 0.7)
```

### Cumulative mass functions

The cumulative mass function of a distribution is the probability that the outcome will be $q$ or less. You can compute the cumulative mass function for the binomial distribution using `pbinom`.



For example, here is the probability of the coin's coming up heads at most once:

```{r}
pbinom(q = 1, size = 10, prob = 0.7)
```

We can check that this works the same way with `dbinom`:

```{r}
dbinom(x = 0, size = 10, prob = 0.7) + dbinom(x = 1, size = 10, prob = 0.7)
```

We indeed got the same number.

We could compute the cumulative mass function for the Bernoulli distribution too -- it would just be kind of boring! (Why?)


### The Guassian distribution

Like we saw in the previous lecture, a Binomial distribution is "bell-shaped" if the number of trials `size` is large:

```{r}
set.seed(0)
exp10 <- rbinom(n = 50000, size = 1000, prob = 0.7)
exp10.df <- data.frame(n.heads = exp10)
ggplot(exp10.df) + geom_bar(mapping = aes(x=n.heads, y = ..prop..)) + scale_x_continuous(limits = c(620, 780))
```

This is the shape of a Gaussian (or Normal) distribution. Specifically, it's (approximately) the shape of a normal distribution with mean $700$ and standard deviation $14.9$. (We'll see how to figure out those numbers later)

Here is how we can generate data from the normal distribution:

```{r}
rnorm(n = 1, mean = 700, sd = 14.9)
rnorm(n = 1, mean = 700, sd = 14.9)
rnorm(n = 1, mean = 700, sd = 14.9)
rnorm(n = 1, mean = 700, sd = 14.9)
```

As you can see, we are now generating real numbers rather than just integers. The Gaussian distribution is a *continuous* distribution -- that means the random process generates real numbers. 

This means that the probability of generating any number in particular is 0 -- there are a *lot* of real numbers. (But we can still compute something that's like a PMF that will tell as how likely a number is to be generated relative to another number -- the function that does that is the probability *density* function, or PDF).

Let's generate (almost) the same histogram as before, but using the Gaussian random number generator:

```{r}
exp.norm <- rnorm(n = 50000, mean = 700, sd = 14.9)
exp.norm.df <- data.frame(r = round(exp.norm))
ggplot(exp.norm.df) + geom_bar(mapping = aes(x=r, y = ..prop..)) +  scale_x_continuous(limits = c(620, 780))
```


We will now use the cumulative probability *density* function for the normal distribution to verify that the probability of getting a number smaller than 700 is 0.5:

```{r}
pnorm(q = 700, mean = 700, sd = 14.9)
```

The same is (approximately) true about the Binomial distribution, which looks very similar:

```{r}
pbinom(q = 700, size = 1000, prob = 0.7)
```

We can verify that the distrbutions are similar by checking the cumulative distribution function values elsewhere:

```{r}
pnorm(q = 680, mean = 700, sd = 14.9)
pbinom(q = 680, size = 1000, prob = 0.7)
```

Finally, we can check that this is actually true in our simulations:

```{r}
mean(exp10 <= 680)
mean(exp.norm <= 680)
```

Pretty close!

### More on the Gaussian distribution

The Guassian distribution (or Bell Curve) is defined by its mean -- the average of the random numbers it generates -- and its standard deviation. The standard deviation is the measure of how "spread out" the numbers are.

The standard deviation of a set of numbers $x$ is
$$sd(x_1, x_2, x_3, ..., x_n) = \sqrt{\frac{ (x_1 - \bar{x})^2 + (x_2 - \bar{x})^2 + ... + (x_n - \bar{x})^2}{n}}.$$

Here, $\bar{x} = \frac{x_1 + x_2 + ... + x_n}{n}$ is the mean of the x's.

The standard deviation is the square root of the average squared distance from a datapoint to the mean. If we just measured the average distance to the mean (instead of the square root of the average squared distance), we'd get a similar quantity that's easier to understand.

In R, we can use `sd` to compute the standard deviation:
```{r}
sd(c(1, 2, 10, 2, 1))    # small
sd(c(100, 2, 10, 2, 1))  # large
```

Let's generate some numbers from a Gaussian distribution:
```{r}
nums <- rnorm(n = 200, mean = 20, sd = 15)
```

We can now estimate the mean and standard deviation:
```{r}
mean(nums)
sd(nums)
```

(Actually, $\sqrt{\frac{n}{n-1}} sd(nums)$ is a better estimate for the standard deviation. We won't get into detail here, but the issue is that we don't assume we know the mean, so we estimate the mean using the mean of the 200 numbers, which won't be exact. So if we just compute the `sd`, we'd be systematically underestimating the true standard deviation).

We can *estimate* many distributions using the Gaussian distribution. You already saw this in action -- we can estimate the Binomial distribution using the Gaussian distribution.


```{r}
set.seed(0)
exp10 <- rbinom(n = 50000, size = 1000, prob = 0.7)
exp10.df <- data.frame(n.heads = exp10)
mean(exp10)
sd(exp10)*sqrt(50000/49999)
```

The Gaussian distribution with mean 700 and standard deviation 14.5 would be a good approximation here, empirically.

### Approximating the Binomial Distrbution using the Guassian distribution

For $size$ (number of tosses) large enough (say over 50) and $prob$ reasonably close to 0.5, we can estimate the binomial distribution using the normal distribution with mean $size\times prob$ and standard deviation $\sqrt{size \times prob \times (1-prob)}$.

This is one application of the Central Limit Theorem (we will not get into this much further.)