02-estimation.Rmd

# Estimation

## Load packages and set plotting theme

```{r, message=FALSE}
library("knitr")
library("kableExtra")
library("MASS")
library("patchwork")
library("tidyverse")
```

```{r}
theme_set(theme_classic() + #set the theme 
            theme(text = element_text(size = 20))) #set the default text size

opts_chunk$set(comment = "",
               fig.show = "hold")
```

## Working with distributions

Every distribution that R handles has four functions. There is a root name, for example, the root name for the normal distribution is `norm`. This root is prefixed by one of the letters here:

```{r, echo=F}
tibble(letter = c("`d`","`p`","`q`","`r`"),
  description = c('for "__density__", the density function (probability function (for _discrete_ variables) or probability density function (for _continuous_ variables))',
                  'for "__probability__", the cumulative distribution function',
                  'for "__quantile__", the inverse cumulative distribution function',
                  'for "__random__", a random variable having the specified distribution'),
  example = c("`dnorm()`", "`pnorm()`", "`qnorm()`", "`rnorm()`")) %>% 
kable() %>% 
kable_styling(bootstrap_options = "striped",
              full_width = F)
```

For the normal distribution, these functions are `dnorm`, `pnorm`, `qnorm`, and `rnorm`. For the binomial distribution, these functions are `dbinom`, `pbinom`, `qbinom`, and `rbinom`. And so forth.

You can get more info about the distributions that come with R via running `help(Distributions)` in your console. If you need a distribution that doesn't already come with R, then take a look [here](https://cran.r-project.org/web/views/Distributions.html) for many more distributions that can be loaded with different R packages. 

### Plotting distributions

Here's an easy way to plot distributions in `ggplot2` using the `stat_function()` function. We take a look at a normal distribution of height (in cm) with `mean = 180` and `sd = 10` (as this is the example we run with in class).

```{r, results = "hold"}
ggplot(data = tibble(height = c(150, 210)),
       mapping = aes(x = height)) +
  stat_function(fun = ~ dnorm(., mean = 180, sd = 10))
```

Note that the data frame I created with `tibble()` only needs to have the minimum and the maximum value of the x-range that we are interested in. Here, I chose `150` and `210` as the minimum and maximum, respectively (which is the mean +/- 3 standard deviations). 

The `stat_function()` is very flexible. We can define our own functions and plot these like here: 

```{r, results='hold'}
# define the breakpoint function 
fun.breakpoint = function(x, breakpoint){
  x[x < breakpoint] = breakpoint
  return(x)
}

# plot the function
ggplot(data = tibble(x = c(-5, 5)),
       mapping = aes(x = x)) +
  stat_function(fun = ~ fun.breakpoint(., breakpoint = 2))
```

Here, I defined a breakpoint function. If the value of `x` is below the breakpoint, `y` equals the value of the breakpoint. If the value of `x` is greater than the breakpoint, then `y` equals `x`. 


### Sampling from distributions

For each distribution, R provides a way of sampling random number from this distribution. For the normal distribution, we can use the `rnorm()` function to take random samples. 

So let's take some random samples and plot a histogram. 

```{r}
# make this example reproducible 
set.seed(1)

# define how many samples to draw 
tmp.nsamples = 100

# make a data frame with the samples
df.plot = tibble(height = rnorm(n = tmp.nsamples, mean = 180, sd = 10))

# plot the samples using a histogram 
ggplot(data = df.plot,
       mapping = aes(x = height)) +
  geom_histogram(binwidth = 1,
                 color = "black",
                 fill = "lightblue") +
  scale_x_continuous(breaks = c(160, 180, 200)) +
  coord_cartesian(xlim = c(150, 210), expand = F)

# remove all variables with tmp in their name 
rm(list = ls() %>% str_subset(pattern = "tmp."))
```

Let's see how many samples it takes to closely approximate the shape of the normal distribution with our histogram of samples. 

```{r}
# make this example reproducible 
set.seed(1)

# play around with this value
# tmp.nsamples = 100
tmp.nsamples = 10000
tmp.binwidth = 1

# make a data frame with the samples
df.plot = tibble(height = rnorm(n = tmp.nsamples,
                                mean = 180,
                                sd = 10)) 

# adjust the density of the normal distribution based on the samples and binwidth 
fun.dnorm = function(x, mean, sd, n, binwidth){
  dnorm(x = x, mean = mean, sd = sd) * n * binwidth
}

# plot the samples using a histogram 
ggplot(data = df.plot,
       mapping = aes(x = height)) +
  geom_histogram(binwidth = tmp.binwidth,
                 color = "black",
                 fill = "lightblue") +
  stat_function(fun = ~ fun.dnorm(x = .,
                                  mean = 180,
                                  sd = 10,
                                  n = tmp.nsamples,
                                  binwidth = tmp.binwidth),
                xlim = c(min(df.plot$height), max(df.plot$height)),
                size = 2) +
  annotate(geom = "text",
           label = str_c("n = ", tmp.nsamples),
           x = -Inf,
           y = Inf,
           hjust = -0.1,
           vjust = 1.1,
           size = 10,
           family = "Courier New") +
  scale_x_continuous(breaks = c(160, 180, 200)) +
  coord_cartesian(xlim = c(150, 210), expand = F)

# remove all variables with tmp in their name 
rm(list = ls() %>% str_subset(pattern = "tmp."))
```

With 10,000 samples, our histogram of samples already closely resembles the theoretical shape of the normal distribution. 

To keep my environment clean, I've named the parameters `tmp.nsamples` and `tmp.binwidth` and then, at the end of the code chunk, I removed all variables from the environment that have "tmp." in their name using the `ls()` function (which prints out all variables in the environment as a vector), and the `str_subset()` function which filters out only those variables that contain the specified pattern.

### Understanding `density()`

First, let's calculate the density for a set of observations and store them in a data frame.

```{r, fig.cap='Density estimation.'}

# calculate density
observations = c(1, 1.2, 1.5, 2, 3)
bandwidth = 0.25 # bandwidth (= sd) of the Gaussian distribution 
tmp.density = density(observations,
        kernel = "gaussian",
        bw = bandwidth,
        n = 512)

# save density as data frame 
df.density = tibble(x = tmp.density$x,
                    y = tmp.density$y) 

df.density %>% 
  head() %>% 
  kable(digits = 3) %>% 
  kable_styling(bootstrap_options = "striped",
                full_width = F)
```

Now, let's plot the density. 

```{r}
ggplot(data = df.density, 
       mapping = aes(x = x, y = y)) +
  geom_line(size = 2) +
  geom_point(data = enframe(observations),
             mapping = aes(x = value, y = 0),
             size = 3)
```

This density shows the sum of the densities of normal distributions that are centered at the observations with the specified bandwidth. 

```{r}
# add densities for the individual normal distributions
for (i in 1:length(observations)){
  df.density[[str_c("observation_",i)]] = dnorm(df.density$x,
                                                mean = observations[i],
                                                sd = bandwidth)
}

# sum densities
df.density = df.density %>%
  mutate(sum_norm = rowSums(select(., contains("observation_"))),
         y = y * length(observations))

df.density %>% 
  head() %>% 
  kable(digits = 3) %>% 
  kable_styling(bootstrap_options = "striped",
                full_width = F)
```

Now, let's plot the individual densities as well as the overall density.

```{r}
# colors of individual Gaussian distributions 
colors = c("blue", "green", "red", "purple", "orange")

# original density 
p = ggplot(data = df.density, aes(x = x, y = y)) +
  geom_line(size = 2)

# individual densities 
for (i in 1:length(observations)){
  p = p + stat_function(fun = "dnorm",
                        args = list(mean = observations[i], sd = bandwidth),
                        color = colors[i])
}

# individual observations 
p = p + geom_point(data = enframe(observations),
             mapping = aes(x = value, y = 0, color = factor(1:5)),
             size = 3,
             show.legend = F) +
  scale_color_manual(values = colors)

# sum of the individual densities
p = p +
  geom_line(data = df.density,
            aes(x = x, y = sum_norm),
            size = 1,
            color = "red",
            linetype = 2)
p # print the figure
```

Here are the same results when specifying a different bandwidth: 

```{r}
# colors of individual Gaussian distributions 
colors = c("blue", "green", "red", "purple", "orange")

# calculate density
observations = c(1, 1.2, 1.5, 2, 3)
bandwidth = 0.5 # bandwidth (= sd) of the Gaussian distribution 
tmp.density = density(observations,
        kernel = "gaussian",
        bw = bandwidth,
        n = 512)

# save density as data frame 
df.density = tibble(
  x = tmp.density$x,
  y = tmp.density$y
) 

# add densities for the individual normal distributions
for (i in 1:length(observations)){
  df.density[[str_c("observation_",i)]] = dnorm(df.density$x,
                                                mean = observations[i],
                                                sd = bandwidth)
}

# sum densities
df.density = df.density %>%
  mutate(sum_norm = rowSums(select(., contains("observation_"))),
         y = y * length(observations))

# original plot 
p = ggplot(data = df.density, aes(x = x, y = y)) +
  geom_line(size = 2) +
  geom_point(data = enframe(observations),
             mapping = aes(x = value,
                           y = 0,
                           color = factor(1:5)),
             size = 3,
             show.legend = F) +
  scale_color_manual(values = colors)

# add individual Gaussians
for (i in 1:length(observations)){
  p = p + stat_function(fun = "dnorm",
                        args = list(mean = observations[i], sd = bandwidth),
                        color = colors[i])
}

# add the sum of Gaussians
p = p +
  geom_line(data = df.density,
            aes(x = x, y = sum_norm),
            size = 1,
            color = "red",
            linetype = 2)
p
```

### Cumulative probability distributions

```{r}
ggplot(data = tibble(height = c(150, 210)),
       mapping = aes(x = height)) +
  stat_function(fun = ~ pnorm(q = ., 
                              mean = 180, 
                              sd = 10)) + 
  labs(x = "height", y = "cumulative probability") + 
  scale_x_continuous(breaks = c(160, 180, 200)) + 
  coord_cartesian(xlim = c(150, 210),
                  ylim = c(0, 1.05),
                  expand = F)
```

Let's find the cumulative probability of a particular value. 

```{r}
tmp.x = 190
tmp.y = pnorm(tmp.x, mean = 180, sd = 10)

print(tmp.y %>% round(3))

# draw the cumulative probability distribution and show the value
ggplot(data = tibble(height = c(150, 210)),
       mapping = aes(x = height)) +
  stat_function(fun = ~ pnorm(q = .,
                              mean = 180, 
                              sd = 10 )) +
  annotate(geom = "point",
           x = tmp.x, 
           y = tmp.y,
           size = 4,
           color = "blue") +
  geom_segment(mapping = aes(x = tmp.x,
                             xend = tmp.x,
                             y = 0,
                             yend = tmp.y),
               size = 1,
               color = "blue") +
  geom_segment(mapping = aes(x = -5,
                             xend = tmp.x,
                             y = tmp.y,
                             yend = tmp.y),
               size = 1,
               color = "blue") +
  labs(x = "height", y = "cumulative probability") + 
  scale_x_continuous(breaks = c(160, 180, 200)) + 
  coord_cartesian(xlim = c(150, 210),
                  ylim = c(0, 1.05),
                  expand = F)

# remove all variables with tmp in their name 
rm(list = str_subset(string = ls(), pattern = "tmp."))
```

Let's illustrate what this would look like using a normal density plot. 

```{r}
ggplot(data = tibble(height = c(150, 210)),
       mapping = aes(x = height)) + 
  stat_function(fun = ~ dnorm(., mean = 180, sd = 10),
                geom = "area",
                fill = "lightblue",
                xlim = c(150, 190)) +
  stat_function(fun = ~ dnorm(., mean = 180, sd = 10),
                size = 1.5) +
  labs(x = "height", y = "density") + 
  scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
  scale_x_continuous(breaks = c(160, 180, 200)) + 
  coord_cartesian(xlim = c(150, 210))
```

### Inverse cumulative distribution

```{r}
ggplot(data = tibble(probability = c(0, 1)),
       mapping = aes(x = probability)) +
  stat_function(fun = ~ qnorm(p = ., 
                              mean = 180,
                              sd = 10)) + 
  labs(y = "height", 
       x = "cumulative probability") + 
  scale_x_continuous(breaks = seq(from = 0, to = 1, by = 0.1)) + 
  scale_y_continuous(limits = c(160, 200)) + 
  coord_cartesian(xlim = c(0, 1.05),
                  expand = F)
```

And let's compute the inverse cumulative probability for a particular value. 

```{r}
tmp.x = 0.3
tmp.y = qnorm(tmp.x, mean = 180, sd = 10)

print(tmp.y %>% round(3))

# draw the cumulative probability distribution and show the value
ggplot(data = tibble(probability = c(0, 1)),
       mapping = aes(x = probability)) +
  stat_function(fun = ~ qnorm(., mean = 180, sd = 10)) +
  annotate(geom = "point",
           x = tmp.x, 
           y = tmp.y,
           size = 4,
           color = "blue") +
  geom_segment(mapping = aes(x = tmp.x,
                             xend = tmp.x,
                             y = 160,
                             yend = tmp.y),
               size = 1,
               color = "blue") +
  geom_segment(mapping = aes(x = 0,
                             xend = tmp.x,
                             y = tmp.y,
                             yend = tmp.y),
               size = 1,
               color = "blue") +
  labs(x = "cumulative probability",
       y = "height") + 
  scale_x_continuous(breaks = seq(from = 0, to = 1, by = 0.1)) + 
  scale_y_continuous(limits = c(160, 200)) + 
  coord_cartesian(xlim = c(0, 1.05),
                  expand = F)

# remove all variables with tmp in their name 
rm(list = str_subset(string = ls(), pattern = "tmp."))
```

### Computing probabilities

#### Via probability distributions

Let's compute the probability of observing a particular value $x$ in a given range. 

```{r}
tmp.lower = 170
tmp.upper = 180

tmp.prob = pnorm(tmp.upper, mean = 180, sd = 10) - 
  pnorm(tmp.lower, mean = 180, sd = 10)

tmp.prob

ggplot(data = tibble(x = c(150, 210)),
       mapping = aes(x = x)) + 
  stat_function(fun = ~ dnorm(., mean = 180, sd = 10),
                geom = "area",
                fill = "lightblue",
                xlim = c(tmp.lower, tmp.upper),
                color = "black",
                linetype = 2) +
  stat_function(fun = ~ dnorm(., mean = 180, sd = 10),
                size = 1.5) +
  labs(x = "height",
       y = "density") + 
  scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
  scale_x_continuous(breaks = c(160, 180, 200)) + 
  coord_cartesian(xlim = c(150, 210))

# remove all variables with tmp in their name 
rm(list = str_subset(string = ls(), pattern = "tmp."))
```

We find that ~34% of the heights are between 170 and 180 cm. 

#### Via sampling

We can also compute the probability of observing certain events using sampling. We first generate samples from the desired probability distribution, and then use these samples to compute our statistic of interest. 

```{r}
# let's compute the probability of observing a value within a certain range 
tmp.lower = 170
tmp.upper = 180

# make example reproducible
set.seed(1)

# generate some samples and store them in a data frame 
tmp.nsamples = 10000

df.samples = tibble(height = rnorm(n = tmp.nsamples, mean = 180, sd = 10))

# compute the probability that sample lies within the range of interest
tmp.prob = df.samples %>% 
  filter(height >= tmp.lower,
         height <= tmp.upper) %>% 
  summarize(prob = n()/tmp.nsamples)

# illustrate the result using a histogram 
ggplot(data = df.samples,
       mapping = aes(x = height)) + 
  geom_histogram(binwidth = 1,
                 color = "black",
                 fill = "lightblue") +
  geom_vline(xintercept = tmp.lower,
             size = 1, 
             color = "red",
             linetype = 2) +
  geom_vline(xintercept = tmp.upper,
             size = 1, 
             color = "red",
             linetype = 2) +
  annotate(geom = "label",
           label = str_c(tmp.prob %>% round(3) * 100, "%"),
           x = 175,
           y = 200,
           hjust = 0.5,
           size = 10) + 
  scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
  labs(x = "height")

# remove all variables with tmp in their name 
rm(list = str_subset(string = ls(), pattern = "tmp."))
```


## Estimating a coin 

```{r, fig.cap='A demonstration of the law of large numbers.'}

# Example taken from here: http://statsthinking21.org/probability.html#empirical-frequency

set.seed(1) # set the seed so that the outcome is consistent
nsamples = 50000 # how many flips do we want to make?

# create some random coin flips using the rbinom() function with
# a true probability of 0.5

df.samples = tibble(trial_number = seq(nsamples), 
                    outcomes = rbinom(nsamples, 1, 0.5)) %>% 
  mutate(mean_probability = cumsum(outcomes) / seq_along(outcomes)) %>% 
  filter(trial_number >= 10) # start with a minimum sample of 10 flips

ggplot(data = df.samples, 
       mapping = aes(x = trial_number, y = mean_probability)) +
  geom_hline(yintercept = 0.5, color = "gray", linetype = "dashed") +
  geom_line() +
  labs(x = "Number of trials",
       y = "Estimated probability of heads") +
  theme_classic() +
  theme(text = element_text(size = 20))
```

## Bayesian inference with the normal distribution

Let's consider the following scenario. You are helping out at a summer camp. This summer, two different groups of kids go to the same summer camp. The chess kids, and the basketball kids. The chess summer camp is not quite as popular as the basketball summer camp (shocking, I know!). In fact, twice as many children have signed up for the basketball camp. 

When signing up for the camp, the children were asked for some demographic information including their height in cm. Unsurprisingly, the basketball players tend to be taller on average than the chess players. In fact, the basketball players' height is approximately normally distributed with a mean of 180cm and a standard deviation of 10cm. For the chess players, the mean height is 170cm with a standard deviation of 8cm. 

At the camp site, a child walks over to you and asks you where their gym is. You gauge that the child is around 175cm tall. Where should you direct the child to? To the basketball gym, or to the chess gym? 

### Analytic solution

```{r}
height = 175

# priors 
prior_basketball = 2/3 
prior_chess = 1/3 

# likelihood  
mean_basketball = 180
sd_basketball = 10

mean_chess = 170
sd_chess = 8

likelihood_basketball = dnorm(height, mean = mean_basketball, sd = sd_basketball)
likelihood_chess = dnorm(height, mean = mean_chess, sd = sd_chess)

# posterior
posterior_basketball = (likelihood_basketball * prior_basketball) / 
  ((likelihood_basketball * prior_basketball) + (likelihood_chess * prior_chess))

posterior_basketball %>% print()
```

### Solution via sampling

Let's do the same thing via sampling. 

```{r}
# number of kids 
tmp.nkids = 10000

# make reproducible 
set.seed(1)

# priors 
prior_basketball = 2/3 
prior_chess = 1/3 

# likelihood functions 
mean_basketball = 180
sd_basketball = 10

mean_chess = 170
sd_chess = 8

# data frame with the kids
df.camp = tibble(kid = 1:tmp.nkids,
                 sport = sample(c("chess", "basketball"),
                                size = tmp.nkids,
                                replace = T,
                                prob = c(prior_chess, prior_basketball))) %>% 
  rowwise() %>% 
  mutate(height = ifelse(test = sport == "chess",
                         yes = rnorm(., mean = mean_chess, sd = sd_chess),
                         no = rnorm(., mean = mean_basketball, sd = sd_basketball))) %>% 
  ungroup

df.camp %>% print()
```

Now we have a data frame with kids whose height was randomly sampled depending on which sport they do. I've used the `sample()` function to assign a sport to each kid first using the `prob = ` argument to make sure that a kid is more likely to be assigned the sport "basketball" than "chess". 

Note that the solution above is not particularly efficient since it uses the `rowwise()` function to make sure that a different random value for height is drawn for each row. Running this code will get slow for large samples. A more efficient solution would be the following: 

```{r}
# number of kids 
tmp.nkids = 100000

# make reproducible 
set.seed(3)

df.camp2 = tibble(
  kid = 1:tmp.nkids,
  sport = sample(c("chess", "basketball"),
                 size = tmp.nkids,
                 replace = T,
                 prob = c(prior_chess, prior_basketball))) %>% 
  arrange(sport) %>% 
  mutate(height = c(rnorm(sum(sport == "basketball"),
                          mean = mean_basketball,
                          sd = sd_basketball),
                    rnorm(sum(sport == "chess"),
                          mean = mean_chess,
                          sd = sd_chess)))
```

In this solution, I take advantage of the fact that `rnorm()` is vectorized. That is, it can produce many random draws in one call. To make this work, I first arrange the data frame, and then draw the correct number of samples from each of the two distributions. This works fast, even if I'm drawing a large number of samples. 

How can we now use these samples to answer our question of interest? Let's see what doesn't work first: 

```{r, eval=F}
tmp.height = 175

df.camp %>% 
  filter(height == tmp.height) %>% 
  count(sport) %>% 
  pivot_wider(names_from = sport, values_from = n) %>% 
  summarize(prob_basketball = basketball/(basketball + chess))
```

The reason this doesn't work is because none of our kids is exactly 175cm tall. Instead, we need to filter kids that are within a certain height range. 

```{r}
tmp.height = 175
tmp.margin = 1

df.camp %>% 
  filter(between(height,
                 left = tmp.height - tmp.margin,
                 right = tmp.height + tmp.margin)) %>% 
  count(sport) %>% 
  pivot_wider(names_from = sport,
              values_from = n) %>% 
  summarize(prob_basketball = basketball/(basketball + chess))
```

Here, I've used the `between()` function which is a shortcut for otherwise writing `x >= left & x <= right`. You can play around with the margin to see how the result changes. 

## Additional resources

### Datacamp

- [Foundations of probability in R](https://www.datacamp.com/courses/foundations-of-probability-in-r)
  
## Session info

Information about this R session including which version of R was used, and what packages were loaded. 

```{r}
sessionInfo()
```