---
title: "Homework: The Idea of Calibration"
tlda-classroom-server: "https://pic.cormorant-matrix.ts.net"
tlda-classroom-assignment: "week1-homework"
execute:
  cache: true
tlda-answer-baseline: "homework-calibration.qmd.support/baseline.txt"
filters:
  - "homework-calibration.qmd.support/answer-placement-warning.lua"
---

```{r}
#| echo: false
required.packages = c("ggplot2", "purrr")
missing.packages = required.packages[!vapply(required.packages, requireNamespace,
                                             quietly = TRUE, FUN.VALUE = logical(1))]
if (length(missing.packages) > 0) {
  install.command = sprintf(
    "install.packages(c(%s))",
    paste(sprintf('"%s"', missing.packages), collapse = ", ")
  )
  stop(
    "This homework needs an R package you do not have installed. Run this in R, then render again:\n\n",
    install.command,
    call. = FALSE
  )
}

library(ggplot2)
library(purrr)
approx = function(x) { round(x, 2) }
theme_set(theme(
  plot.background = element_rect(fill = "transparent", colour = NA),
  panel.background = element_rect(fill = "transparent", colour = NA),
  legend.background = element_rect(fill = "transparent", colour = NA),
  legend.box.background = element_rect(fill = "transparent", colour = NA),
  legend.key = element_rect(fill = "transparent", colour = NA)
))
hawks.six = data.frame(
  roll   = 1:6,
  player = c("Dejounte Murray", "Trae Young", "John Collins",
             "Saddiq Bey", "De'Andre Hunter", "Onyeka Okongwu"),
  points = c(1515, 1914, 931, 1062, 1029, 791)
)
```

::: {.callout-note title="On the timing of this one"}
The first two sections are meant to be done *before* Thursday's class --- they are what that
class builds on. The rest is easier once you have been to it. The due date is Thursday at
midnight, which does not really give you time to do everything after the class, so do the
early parts first and don't save them.
:::

This homework is about calibrating interval estimates. We will cheat by working with
populations we know. That lets us build the sampling distribution, choose an interval width,
and check how often the resulting intervals cover the truth.

We start with a die, where we can do every calculation by hand. Then we use the same
procedure on NBA data.

# A Die

Roll a six-sided die. If you don't have one, ask me and I'll give you one.

::: {#exr-die-histogram .callout-exercise}
Roll it thirty times and write down how often each face came up.

Then, **by hand** on paper, draw the histogram of those thirty rolls --- one bar per face,
its height the number of times that face came up. Photograph it and put it in your answer
along with your six counts.

Each of the six faces is equally likely --- the probability of rolling any particular one is
$1/6$. Does that mean the six bars should come out the same height?
:::


::: {#ans-exr-die-histogram .callout-answer .callout-note title="Template"}

*(your six counts, your photo, and your answer)*

:::

## A population of six sides

Think of the sides of the die as a population. Write $j=1,\ldots,6$ for a side's position
and $y_j$ for the value written on it. For this die, $y_j=j$.

```{r}
#| echo: false
knitr::kable(data.frame(j = 1:6, y = 1:6),
             col.names = c("side ($j$)", "value ($y_j$)"))
```

The first roll, $J_1$, chooses a member of that population, every one equally likely. The
value we observe is $Y_1=y_{J_1}$.

probability | value of $Y_1$
----|---
1/6 | 1
1/6 | 2
1/6 | 3
1/6 | 4
1/6 | 5
1/6 | 6

This table is the *distribution* of $Y_1$.

Suppose we want to know the probability that a roll is 3 or more. Nobody looks
this up in a book. They count. Four of the six outcomes are 3 or more, so the
probability is 4/6.

The probability it doesn't happen---the *complementary probability*---is 2/6.

probability | value
----|---
2/6 | 0
4/6 | 1

This is the distribution of the indicator $1_{\ge 3}(Y_1)$, which is defined like
this:

$$
1_{\ge 3}(y) = \begin{cases} 1 & \text{ if } y \ge 3 \\ 0 & \text{ otherwise.} \end{cases}
$$

Counting worked because two different faces cannot come up on the same roll. The probability
of getting one face *or* another is therefore the sum of their probabilities. That is
*disjoint additivity* of probabilities. We'll talk about this later in the semester when we
get more formal about probability.

## Marginalization

We just worked this out by counting. Now we're going to formalize that idea in a way that
generalizes pretty well.

Instead of going directly from the probability distribution of one thing to the probability
distribution of another, we talk about the probability distribution of both things together
--- the probability distribution of the pair.

For a six-valued die and a two-valued indicator, there are $6 \times 2 = 12$ possible pairs.
Thankfully, this is an easy case. Our indicator $1_{\ge 3}(Y_1)$ is a *function* of the die
value $Y_1$. Whenever we know the value of $Y_1$, we also know the value of
$1_{\ge 3}(Y_1)$. If we wanted to put the other six pairs into our table, we could, but we
know that they don't happen. Their probability is zero. So practically speaking, we can add
a function to our probability table for free. We just add a column with the function's value
for each die value.

probability | $Y_1$ | $1_{\ge 3}(Y_1)$
--- | --- | ---
1/6 | 1 | 0
1/6 | 2 | 0
1/6 | 3 | 1
1/6 | 4 | 1
1/6 | 5 | 1
1/6 | 6 | 1

Then we get the distribution of the indicator by summing the probabilities of the rows with
the same value in the indicator column. This is called *marginalization*.

## Estimation

We ask a question about the six sides of the die: what fraction are 3 or more? That is, if
we call the sides $y_1,\ldots,y_6$, we want this:

$$
\underset{\text{estimand}}{\theta} = \frac{1}{6}\sum_{j=1}^6 1_{\ge 3}(y_j) = \frac{4}{6}.
$$

We're going to estimate it by rolling $n$ times and calling the values $Y_1,\ldots,Y_n$.
This is a sample drawn **with
replacement**: the same value can appear twice because the die does not know what it did last
time. The fraction of rolls that are 3 or more is our *estimate*:

$$
\underset{\text{estimate}}{\hat\theta_n} = \frac{1}{n}\sum_{i=1}^n 1_{\ge 3}(Y_i).
$$

Here we're following a convention in which we write our targets as some greek letter. Often
$\theta$. And add a 'hat', the thing above the letter theta in $\hat\theta$, to indicate an
estimator of that thing.

When we repeat the sample, $\hat\theta_n$ changes. Its *sampling distribution* lists the values
it can take and how often each one occurs. We can work it out exactly here because the
population is small enough to count.

### Sample Size One

::: {#exr-count-n1 .callout-exercise}
Use your thirty rolls to approximate the sampling distribution of $\hat\theta_1$. Compare
your approximation with the exact distribution in the indicator table above.
:::


::: {#ans-exr-count-n1 .callout-answer .callout-note title="Template"}

*(your approximation and comparison)*

:::

### Sample Size Two

::: {#exr-count-n2 .callout-exercise}
Now take $n=2$. Make a table of the possible pairs of die sides and the value of
$\hat\theta_2$ for each pair, then marginalize the table to get the sampling distribution
of $\hat\theta_2$. Write out the calculation you used.

Could you have made this easier on yourself? *Hint*. What did you learn about addition in
middle school?
:::


::: {#ans-exr-count-n2 .callout-answer .callout-note title="Template"}

*(your table and calculation)*

:::

::: {#exr-compare-sample-sizes .callout-exercise}
The opening exercise gave you thirty values of $\hat\theta_1$. Now roll two dice thirty times
and calculate $\hat\theta_2$ each time. Compare the two distributions. What changed when the
sample size went from one to two?
:::


::: {#ans-exr-compare-sample-sizes .callout-answer .callout-note title="Template"}

*(your thirty values and your comparison)*

:::

::: {#exr-width-by-hand .callout-exercise}
Use your thirty values of $\hat\theta_2$ to draw thirty interval estimates, each centred at
one value. Give them all the same width $w$, so the interval centred at
$\hat\theta_2$ runs from $\hat\theta_2-w/2$ to $\hat\theta_2+w/2$. Choose the smallest $w$
for which at least twenty of the thirty intervals contain $\theta=2/3$.

Start at $\theta$ and count outward through your thirty values. Draw the intervals, report
their width $w$, and count their coverage.
:::


::: {#ans-exr-width-by-hand .callout-answer .callout-note title="Template"}

*(your thirty intervals, width, and coverage)*

:::

This function does what you just did by hand. It calculates every draw's distance from the
centre, puts those distances in order, and takes the first one that reaches the target
coverage. Then it doubles that distance to turn a half-width into a width.

```{r}
#| label: width-code
#| echo: true
width = function(draws, center = mean(draws), alpha = .05) {
    2 * quantile(abs(draws - center), 1 - alpha, type = 1)
}
```

# Simulating

Counting worked because there were 36 pairs. At $n=3$ there are 216, at $n=5$ there are
7776, and the assignment is not going to ask you to count those.

We are not going to ask you to do these larger calculations by hand. The code does the same
things you just did: build the sampling distribution, then count outward from its centre until
the interval reaches its target coverage.

::: {#exr-simulate-n2 .callout-exercise}
You worked out the sampling distributions and calibrated the interval estimates by hand.
Now use the code below to check your work. We are only asking you to fill in the estimator
and the estimand.

For $n=2$, calibrate the interval width and check the coverage. Compare them with
the interval width and coverage you calculated by hand.
:::


::: {#ans-exr-simulate-n2 .callout-answer .callout-note title="Template"}

```{r}
#| eval: false
#| echo: true
population = 1:6
cutoff = 3

estimator = function(sample, cutoff) {
  NA                              # <1>
}

theta = NA                        # <2>

populations = rep(list(population), 10000)

n = 2
theta.hats = populations |> map_vec(function(population) {
  sample = sample(population, n, replace = TRUE)
  estimator(sample, cutoff)
})

ggplot() + geom_bar(aes(x = theta.hats, y = after_stat(prop))) +
  labs(x = expression(hat(theta)), y = "")

w = width(theta.hats, center = theta, alpha = 1/3)
w
mean(theta.hats - w/2 <= theta & theta <= theta.hats + w/2)
```
1. Fill in the estimator: the fraction of the sample at or above the cutoff.
2. Calculate the estimand for this population and cutoff.

The rest of the block simulates and plots the sampling distribution, calibrates the interval
width, and checks the resulting coverage.

:::

::: {#exr-simulate-n .callout-exercise}
Now run it at $n=10$ and $n=50$, and put the histograms for $n=2$, $n=10$, and $n=50$
next to each other.

What happens to the sampling distribution as $n$ grows? Say what happens to where it sits
and what happens to how wide it is, and say which of those two you would have predicted.
:::


::: {#ans-exr-simulate-n .callout-answer .callout-note title="Template"}

```{r}
#| echo: true
# your three simulations and your three histograms
```

*(your answer here)*

:::

# NBA Data

## A population of six players

Now do the same thing with a different population. Here are the Atlanta Hawks' starters and
their sixth man, along with the points each scored over the 2023 season.

```{r}
#| echo: false
knitr::kable(hawks.six, col.names = c("roll ($j$)", "player", "points ($y_j$)"))
```

Rolling the die chooses one player, every one equally likely. Now $Y_1$ is the chosen
player's point total.

::: {#exr-die-estimate .callout-exercise}
Add an indicator column for whether the player scored at least 1000 points. Then marginalize
the six-row table to find the distribution of that indicator.
:::


::: {#ans-exr-die-estimate .callout-answer .callout-note title="Template"}

*(your table)*

:::

::: {#exr-hawks-n2 .callout-exercise}
Now work out the sampling distribution of $\hat\theta_2$, where the indicator is 1 when the
selected player scored at least 1000 points. Write out the calculation you used.
:::


::: {#ans-exr-hawks-n2 .callout-answer .callout-note title="Template"}

*(your table and calculation)*

:::

You have now repeated the hand calculation with the Hawks' six. Now repeat the code check
using the estimator you defined for the die. The sampling, replication, plot, interval width,
and coverage calculation are supplied.

::: {#exr-hawks-simulation .callout-exercise}
Use the code below to simulate the sampling distribution of $\hat\theta_2$ for the Hawks'
six. Compare it with the distribution you calculated by hand. Then use the simulated
distribution to calibrate the interval width and check its coverage.
:::


::: {#ans-exr-hawks-simulation .callout-answer .callout-note title="Template"}

```{r}
#| eval: false
#| echo: true
population = hawks.six$points
cutoff = 1000
theta = estimator(population, cutoff)

populations = rep(list(population), 10000)
n = 2
theta.hats = populations |> map_vec(function(population) {
  sample = sample(population, n, replace = TRUE)
  estimator(sample, cutoff)
})

ggplot() + geom_bar(aes(x = theta.hats, y = after_stat(prop))) +
  labs(x = expression(hat(theta)), y = "")

w = width(theta.hats, center = theta, alpha = 1/3)
w
mean(theta.hats - w/2 <= theta & theta <= theta.hats + w/2)
```
This calls the estimator you defined above with the Hawks' population and a cutoff of 1000.

:::

## Two Samples

A mutual friends gives us some data describing the performance of NBA players in the 2023 season. Actually, they give us different data---they give me 'sample 1' and you 'sample 2'. They say these are samples of size 100 drawn with replacement from the population of all players that played in the league that year.  

```{r}
#| echo: true
my.sample = read.csv("https://qtm285-1.github.io/assets/data/nba_sample_1.csv")
your.sample = read.csv("https://qtm285-1.github.io/assets/data/nba_sample_2.csv")
```

Theres a fair amount of information about the players in our samples. Here's a list of the variables we have.

- Player - The player's name.
- Team: The team whose this player is playing for this season in abbreviated term.
- Age:  The Age of the player.
- GP: The number of games that the player has played in this season.
- W: The number of games won in which the player has played.
- L: The numer of games lost in which the player has played.
- Min: The minutes the player has played for this season.
- PTS: The number of points made by the player.

I'll work through my sample first. Then you'll repeat the same analysis with yours.

## An Estimation Problem

We ask the same question we asked about the Hawks' six, on the whole league: what fraction
of players scored at least 1000 points?

My estimate, for what it's worth, is `r approx(mean(my.sample$PTS >= 1000))`.
```{r}
#| echo: true
Y = my.sample$PTS
n = length(Y)

theta.hat = mean(Y >= 1000)
theta.hat
```

To put it in context, I plot my points around the thousand-point line.

```{r}
#| echo: true
#| fig-alt: "Points scored by each player in sample 1, with a red horizontal line at one thousand points."
library(ggplot2)

scale_y = scale_y_continuous(breaks=seq(0,2000,by=500), limits=c(0,2300))
scale_x = scale_x_continuous(breaks=seq(0,2000,by=500), limits=c(0,2300))
no_labs = labs(x='', y='')

points.plot = ggplot(data.frame(i = 1:n, y = Y)) +
    annotate("rect", xmin = -Inf, xmax = Inf, ymin = 1000, ymax = Inf,
             alpha = .1, fill = "red") +
    geom_point(aes(x = i, y = y), alpha = .4) +
    geom_hline(aes(yintercept = 1000), color = "red") +
    scale_y + no_labs 

points.histogram = ggplot(data.frame(y = Y)) +
    annotate("rect", xmin = 1000, xmax = Inf, ymin = -Inf, ymax = Inf,
             alpha = .1, fill = "red") +
    geom_histogram(aes(x = Y, y = after_stat(density)), bins = 20, alpha = .2) +
    geom_vline(aes(xintercept = 1000), color = "red") +
    scale_x + no_labs

points.plot
```

Here the x-coordinate is irrelevant. I've plotted pairs $(i,Y_i)$ where $i=1 \ldots n$ counts out the observations in our sample and $Y_i$ is the number of points the $i$th player in my sample scored. The red horizontal line is the thousand-point mark, so my estimate $\hat\theta$ is the fraction of the dots sitting above it.

To make it a little easier to see the distribution of season point totals in my sample, I plot a histogram and annotate it with the same line---still red, but now vertical to fit into the layout of the histogram. The shaded region is where the indicator equals 1. The area of the histogram bars inside that region approximates $\hat\theta$; when a bin crosses the thousand-point line, the histogram has grouped together points on both sides.

```{r}
#| echo: true
#| fig-alt: "Histogram of points scored by players in sample 1, with a red vertical line at one thousand points."
points.histogram = ggplot(data.frame(y = Y)) +
    annotate("rect", xmin = 1000, xmax = Inf, ymin = -Inf, ymax = Inf,
             alpha = .1, fill = "red") +
    geom_histogram(aes(x = Y, y=after_stat(density)), bins=20, alpha = .2) +
    geom_vline(aes(xintercept = 1000), color = "red") +
    scale_x + no_labs
points.histogram
```

To make these fit together nicely, I rotate my histogram and lay my two plots out side by side so that the bins of the histogram align with the individual players' points.^[When you do this, it's important to make sure the grid lines in your two plots line up. To do that, be sure to set the same limits on the x-axis in both plots and make sure axis labels aren't shifting things up or down. You can do the former by setting `breaks` and `limits` in `scale_*_continuous` and the latter by saying `labs(x='', y='')` to turn axis-labels off.] 


```{r}
#| echo: true
#| layout-ncol: 2
#| fig-show: hold
#| fig-width: 5
#| fig-alt:
#| - "Rotated histogram of sample 1's points, with a red line at one thousand."
#| - "Points scored by each player in sample 1, aligned with the rotated histogram, with a red line at one thousand."

points.histogram + coord_flip()                     # <1>
points.plot + theme(axis.text.y = element_blank())  # <2>
```
1. This `coord_flip()` call rotates the histogram plot so the bins are horizontal.
2. The `theme(...)` call turns off the vertical 'tick labels' 0, 500, ... in the plot on the right. They're redundant because our grids are aligned and we've already got the on the left.

::: {#exr-your-estimate .callout-exercise}

Replicate my analyis using your sample. Report your estimate as well as your two plots.

*Hint*. You should be able to do this with a one-line change to my code. 
:::


::: {#ans-exr-your-estimate .callout-answer .callout-note title="Template"}

```{r}
#| echo: true
# your one-line change and analysis go here
```

:::


My estimate and yours are fairly different, even though our friend says both samples came
from the same population. That difference could just be sampling variation. To find out, we
ask our friend for the data on every player in the league.

Knowing the whole population lets us simulate what estimates from samples of 100 usually
look like. Then we can see whether your estimate is a plausible result of sampling from this
population.

```{r}
#| echo: true

our.pop = read.csv("https://qtm285-1.github.io/assets/data/nba_population.csv") 
```

::: {#exr-sampling-distribution .callout-exercise}
To get a sense of the sampling distribution of the estimator you and I are using, draw 10,000 samples of size $n=100$ with replacement from this population, and get, for each sample, the fraction of the players in it who scored at least 1000 points. That approximates our sampling distribution. Plot the fractions you get in a histogram to get a visual sense and then answer this: how weird is your estimate? Do you think it's likely that your sample was really drawn with replacement from this population?
:::



::: {#ans-exr-sampling-distribution .callout-answer .callout-note title="Template"}

Here is the same interface you used for the die. This time, apply the estimator you already
defined to the league's points data. The rest of the block draws the samples, collects the
estimates, and plots the sampling distribution.

```{r}
#| eval: false
#| echo: true

population = our.pop$PTS
cutoff = 1000
theta = estimator(population, cutoff)

populations = rep(list(population), 10000)
estimates = populations |> map_vec(function(population) {
  n = 100
  sample = sample(population, n, replace = TRUE)
  estimator(sample, cutoff)
})

sampling.distribution.histogram = ggplot() +
  geom_histogram(aes(x = estimates, y = after_stat(density)), bins = 20, alpha = .2)

sampling.distribution.histogram +
  geom_point(aes(x = mean(your.sample$PTS >= 1000), y = .003), color = "red")
```

:::

# Two Populations

You have now calculated the fraction of players scoring at least 1,000 points in two
populations: the league and the Hawks’ six players.

One thing to know while you look at them: the Hawks finished that season 41-41.

::: {#exr-two-populations .callout-exercise}
Compare what you found for the Hawks' six with what you found for the league. Why do you
think they are different? Is there another comparison you should be making? If so, go for
it. You have the code.
:::


::: {#ans-exr-two-populations .callout-answer .callout-note title="Template"}

```{r}
#| echo: true
# your comparison goes here
```

*(your answer here)*

:::
