probsim_tv.Rmd

---
title: "Regression and Other Stories: Simulation"
author: "Andrew Gelman, Jennifer Hill, Aki Vehtari"
date: "`r Sys.Date()`"
output:
  github_document:
    toc: true
---
Tidyverse version by Bill Behrman.

Simulation of probability models. See Chapter 5 in
Regression and Other Stories.

-------------

```{r, message=FALSE}
# Packages
library(tidyverse)

# Parameters
  # Common code
file_common <- here::here("_common.R")

#===============================================================================

# Run common code
source(file_common)
```

# 5 Simulation

## 5.1 Simulation of dicrete probability models

### How many girls in 400 births?

```{r}
n_births <- 400
n_sims <- 1e4

prob_girl <- 0.488
```

One observation of `r n_births` single births.

```{r}
set.seed(151)

rbinom(n = 1, size = n_births, prob = prob_girl)
```

`r format(n_sims, big.mark = ",")` observations of `r n_births` single births.

```{r}
set.seed(151)

sim_girl_1 <- 
  tibble(n_girls = rbinom(n = n_sims, size = n_births, prob = prob_girl))
```

Number of girls born in `r format(n_sims, big.mark = ",")` observations of `r n_births` single births.

```{r}
sim_girl_1 %>% 
  ggplot(aes(n_girls)) +
  geom_histogram(binwidth = 2) +
  labs(
    title = 
      str_glue(
        "Number of girls born in {format(n_sims, big.mark = ',')} observations of {n_births} single births"
      ),
    x = "Number of girls born",
    y = "Count"
  )
```

`r format(n_sims, big.mark = ",")` observations of `r n_births` single and twin births.

```{r}
set.seed(151)

prob_fraternal <- 1 / 125
prob_identical <- 1 / 300
prob_single <- 1 - prob_fraternal - prob_identical
prob_single_girl <- 0.488
prob_fraternal_girl <- 0.495
prob_identical_girl <- 0.495

n_girls <- function(n_single, n_fraternal, n_identical) {
  n_single_girls <- rbinom(n = 1, size = n_single, prob = prob_single_girl)
  v <- 
    rmultinom(
      n = 1,
      size = n_fraternal,
      prob = 
        c(
          (1 - prob_fraternal_girl)^2,
          2 * (1 - prob_fraternal_girl) * prob_fraternal_girl,
          prob_fraternal_girl^2
        )
    ) %>% 
    t() %>% 
    as_tibble(.name_repair = ~ c("n_girls_0", "n_girls_1", "n_girls_2"))
  n_fraternal_girls <- v$n_girls_1 + 2L * v$n_girls_2
  n_identical_girls <- 
    2L * rbinom(n = 1, size = n_identical, prob = prob_identical_girl)
  n_single_girls + n_fraternal_girls + n_identical_girls
}

sim_girls_2 <- 
  tibble(
    n_girls =
      rmultinom(
        n = n_sims,
        size = n_births,
        prob = c(prob_single, prob_fraternal, prob_identical)
      ) %>% 
      t() %>% 
      as_tibble(
        .name_repair = ~ c("n_single", "n_fraternal", "n_identical")
      ) %>% 
      pmap_int(n_girls)
  )
```

Number of girls born in `r format(n_sims, big.mark = ",")` observations of `r n_births` single and twin births.

```{r}
sim_girls_2 %>% 
  ggplot(aes(n_girls)) +
  geom_histogram(binwidth = 2) +
   labs(
    title = 
      str_glue(
        "Number of girls born in {format(n_sims, big.mark = ',')} observations of {n_births} single and twin births"
      ),
    x = "Number of girls born",
    y = "Count"
  )
```

## 5.2 Simulation of continuous and mixed discrete / continuous models

```{r}
mean <- 3
sd <- 0.5
```

Normal distribution with mean `r mean` and standard deviation `r sd`.

```{r}
set.seed(660)

v <- tibble(x = rnorm(n_sims, mean = mean, sd = sd))

v %>% 
  ggplot(aes(x)) +
  geom_histogram(binwidth = 0.1) +
  labs(
    title =
      str_glue(
        "Normal distribution with mean {mean} and standard deviation {sd}"
      )
  )
```

Corresponding log-normal distribution.

```{r}
v %>% 
  ggplot(aes(exp(x))) +
  geom_histogram(binwidth = 2) +
  labs(title = "Corresponding log-normal distribution")
```

```{r}
size <- 20
prob <- 0.6
```

Binomial distribution with `r size` trials and probability `r prob`.

```{r}
set.seed(136)

v <- tibble(x = rbinom(n_sims, size = size, prob = prob))

v %>% 
  ggplot(aes(x)) + 
  geom_bar() +
  scale_x_continuous(breaks = seq(0, size), minor_breaks = NULL) +
  coord_cartesian(xlim = c(0, size)) +
  labs(
    title =
      str_glue(
        "Binomial distribution with {size} trials and probability {prob}"
      )
  )
```

```{r}
lambda <- 5
```

Poisson distribution with mean `r lambda`.

```{r}
set.seed(185)

v <- tibble(x = rpois(n_sims, lambda = lambda))

v %>% 
  ggplot(aes(x)) +
  geom_bar() +
  scale_x_continuous(breaks = seq(0, max(v$x)), minor_breaks = NULL) +
  labs(title = str_glue("Poisson distribution with mean {lambda}"))
```

Simulation parameters on probabilities and heights of adults.

```{r}
n_sims <- 1000

prob_men <- 0.48
height_men_mean <- 69.1
height_men_sd <- 2.9
height_women_mean <- 64.5
height_women_sd <- 2.7
```

Mean height of 10 randomly chosen adults.

```{r}
height_mean <- function(n = 10) {
  if_else(
    rbinom(n, size = 1, prob = prob_men) == 1,
    rnorm(n, mean = height_men_mean, sd = height_men_sd),
    rnorm(n, mean = height_women_mean, sd = height_women_sd)
  ) %>% 
    mean()
}
```

One observation.

```{r}
set.seed(815)

height_mean()
```

`r format(n_sims, big.mark = ",")` observations.

```{r}
set.seed(815)

sim_height_mean <- tibble(height_mean = map_dbl(seq_len(n_sims), height_mean))
```

Distribution of the mean height of 10 adults.

```{r}
sim_height_mean %>%
  ggplot(aes(height_mean)) +
  geom_histogram(binwidth = 0.1) +
  labs(title = "Distribution of the mean height of 10 adults")
```

Maximum height of 10 randomly chosen adults.

```{r}
height_max <- function(n = 10) {
  if_else(
    rbinom(n, size = 1, prob = prob_men) == 1,
    rnorm(n, mean = height_men_mean, sd = height_men_sd),
    rnorm(n, mean = height_women_mean, sd = height_women_sd)
  ) %>% 
    max()
}
```

One observation.

```{r}
set.seed(815)

height_max()
```

`r format(n_sims, big.mark = ",")` observations.

```{r}
set.seed(815)

sim_height_max <- tibble(height_max = map_dbl(seq_len(n_sims), height_max))
```

Distribution of the maximum height of 10 adults.

```{r}
sim_height_max %>%
  ggplot(aes(height_max)) +
  geom_histogram(binwidth = 0.5) +
  scale_x_continuous(breaks = scales::breaks_width(1)) +
  labs(title = "Distribution of the maximum height of 10 adults")
```