---
title: "Homework: RVs, Markov's Inequality, Chebyshev's Inequality, Biased Estimators"
format:
  html:
    theme: [default, "class.scss"]
    callout-appearance: simple
    callout-icon: false
tlda-answer-baseline: "homework-markov-chebyshev.qmd.support/baseline.txt"
filters:
  - "solution-callout.lua"
  - "homework-markov-chebyshev.qmd.support/answer-placement-warning.lua"
---

<!-- Sources:
  Section 1 (Markov's Inequality and Consistency): homework/enrichment-probability.qmd lines 39-81
  Section 2 (Proving Markov's Inequality): homework/enrichment-probability.qmd lines 58-69
  Section 3 (Markov's Inequality and Interval Estimation): homework/enrichment-probability.qmd lines 83-177
  Section 4 (Chebyshev's Inequality): Already derived in Section 3 (applying Markov to the square). Just needs naming.
  Section 5 (Coverage with NBA data): homework/enrichment-probability.qmd lines 131-355 — pull here
  Section 6 (Bias-variance tradeoff): homework/homework2.qmd lines 40-46 (decomposition) + lines 75-110 (visualizations)
  Section 7 (Biased Estimators — Introduction): homework/week2b-homework.qmd lines 19-25
  Section 8 (Using Prior Information): homework/week2b-homework.qmd lines 28-98
  Section 9 (Visualizing the Impact of Prior Information): homework/week2b-homework.qmd lines 101-177
  GAPs: 2 (simulation exercises)
-->

{{< include shared-code.qmd >}}
{{< include ga-turnout.qmd >}}

::: hidden
$$
\DeclareMathOperator{\E}{E}
\DeclareMathOperator{\Var}{V}
\DeclareMathOperator{\P}{P}
\DeclareMathOperator{\sd}{sd}
\DeclareMathOperator{\bias}{bias}
\newcommand{\thetaprior}{\theta^{\text{prior}}}
\newcommand{\nprior}{n^{\text{prior}}}
\newcommand{\yprior}{y^{\text{prior}}}
$$
:::


# Markov's Inequality and Consistency

Convergence in mean square *implies* convergence in probability. Let's use *Markov's inequality* to see why.

## Markov's Inequality

Markov's inequality says that for a non-negative random variable $X$ and any $t > 0$,
$$
P(X \ge t) \le \frac{\E[X]}{t}.
$$

The usual proof of Markov's inequality is based on a few simple observations.

1. The expectation of the indicator variable $1_{\ge t}(X)$ is the probability that $X$ exceeds $t$, $P(X \ge t)$.
2. If we have some function of $u$ that's always larger than $1_{\ge t}$, i.e.,
   one satisfying $u_t(x) \ge 1_{\ge t}(x)$ for all $x$, we know that $\E[u_t(X)] \ge \E 1_{\ge t}(X)$
   for any random variable $X$. If it's always larger for non-negative $x$, then $\E[u_t(X)] \ge \E[1_{\ge t}(X)]$ for any non-negative random variable $X$.
3. The function $u_t(x)=x/t$ is such a function.^[If you're not convinced, sketch the two functions on the same axes. Sketching usually helps.]

::: {.callout-exercise}
## Proving Markov's Inequality
Prove Markov's inequality.
:::

::: {#ans-q1 .callout-answer .callout-note title="Template"}

*(your answer here)*

:::

Often, instead of using this to bound the random variable we're interested in directly,
e.g. $X=|\hat\theta - \theta|$, we use it to bound the random variable's square.
$|X| \ge \epsilon$ if and only if $X^2 \ge \epsilon^2$, so the probability that $|X| \ge \epsilon$ is the same as the probability that $X^2 \ge \epsilon^2$. Applying Markov's inequality to the random variable $X^2$ gives us a bound in terms of $X$'s mean square which, in the specific case
that $X$ is $|\hat\theta-\theta|$, is the *mean squared error* of the estimator $\hat\theta$, $\RMSE(\hat\theta)^2=\E[(\hat\theta-\theta)^2]$

$$
P(X \ge \epsilon) = P(X^2 \ge \epsilon^2) \leq \frac{\E[X^2]}{\epsilon^2} \qqtext{ e.g.}
P(|\hat\theta - \theta| \ge \epsilon) = P((\hat\theta - \theta)^2 \ge \epsilon^2) \leq \frac{\E[(\hat\theta - \theta)^2]}{\epsilon^2}
$$

This tells us that, if the root-mean-squared error of $\hat\theta$ goes to zero, then the probability that $\hat\theta$ is any distance $\epsilon$ away from $\theta$ goes to zero, i.e., consistency in mean-square implies consistency in probability.

::: {#exr-markov-simulation .callout-exercise}
### Markov in Action

Draw $n=10{,}000$ observations from an $\text{Exponential}(1)$ distribution (use `rexp(10000)`).

a. For $\epsilon = 2, 3, 5$, compute the fraction of observations that exceed $\epsilon$. These are estimates of $P(X \ge \epsilon)$.
b. For each $\epsilon$, compute the Markov bound $\E[X]/\epsilon = 1/\epsilon$.
c. Compare your empirical tail probabilities to the Markov bounds. How conservative is the bound?
d. Now try the same thing with $X^2$ in place of $X$. That is, square your observations and compare the fraction exceeding $\epsilon^2$ to $\E[X^2]/\epsilon^2$. Is this version of the bound tighter or looser?
:::

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

*(your answer here)*

:::


# Markov's Inequality and Interval Estimation

So far, when we've calibrated interval estimates using our estimator's standard deviation, we've relied on normal approximation. In effect, we've been using a formula for $P(\lvert\hat\theta - \theta\rvert \le \epsilon)$ that's accurate when $\hat\theta$ has a normal distribution and close enough when its distribution is close enough to normal. In this problem, we're going to think about doing without this reliance on approximate normality.

Let's consider $\hat\theta$, an unbiased estimator of $\theta$ with standard deviation $\sigma$, so the
normal approximation to the distribution of $\hat\theta-\theta$ has the density $f_{0,\sigma}(x)$ below.

$$
P\p*{|\hat\theta - \theta| \le \epsilon}
\approx \int_{-\epsilon}^{\epsilon} f_{0,\sigma}(x) dx
\qfor f_{0, \sigma}(x) = \frac{1}{\sqrt{2\pi}\sigma} e^{-x^2/2\sigma^2}
$$

The reason we've been talking about interval estimators of the form $\hat\theta \pm 1.96 \sigma$ is that, if this approximation were perfect,
these interval estimators would have 95% coverage. That is, it'd be true that $P(|\hat\theta-\theta| \le 1.96 \sigma) = .95$.
And if the approximation is pretty good, we should still expect coverage close to that. But suppose we're not confident that it is.
Markov's inequality allows us to calibrate interval estimates in terms of our estimator's standard deviation
*without any caveats* about its sampling distribution being approximately normal. Let's give it a shot.

::: {#exr-markov-interval  .callout-exercise}
Let $\hat\theta$ be an unbiased estimator of $\theta$ with standard deviation $\sigma$.
By applying Markov's inequality to $|\hat\theta-\theta|^2$, find a *lower bound* on $P(|\hat\theta-\theta| \le t\sigma)$ as a function of $t$.
For what choice of $t$ is this bound equal to $.95$?  Explain how you can use this to get an interval estimate $\hat\theta \pm t\sigma$
with a coverage probability of at least 95%.
:::

::: {.callout-tip collapse="true"}
### From Upper Bounds to Lower Bounds.
If you have an upper bound on $P(X > \epsilon)$, then you have a lower bound on $P(X \le \epsilon)$.
$P(X \le \epsilon)=1-P(X > \epsilon)$ because $X > \epsilon$ if and only if it is not the case that $X \le \epsilon$,
so if $P(X > \epsilon) \le u$, then $P(X \le \epsilon) = 1-P(X > \epsilon) \ge 1-u$. And if you're troubled that
Markov's inequality says $\ge$ instead of $>$, don't be. You have a bound on $P(X > \epsilon)$ too.
$P(X > \epsilon) \le  P(X\ge \epsilon)$ because $X \ge \epsilon$ whenever $X > \epsilon$.
:::

::: {#ans-exr-markov-interval .callout-answer .callout-note title="Template"}

*(your answer here)*

:::

::: {#exr-markov-vs-normal-intervals .callout-exercise}
### Markov vs Normal Intervals

Consider estimating the mean of an $\text{Exponential}(1)$ distribution using the sample mean $\bar X$ from a sample of size $n=25$.

a. Run 10,000 simulations. In each, draw a sample of size 25, compute $\bar X$ and $s = \text{sd}(X)/\sqrt{n}$, and construct two 95% intervals:
   - **Normal interval**: $\bar X \pm 1.96s$
   - **Markov interval**: $\bar X \pm \sqrt{20}\, s$
b. What fraction of each type of interval contains the true mean $\theta=1$? (This is the empirical coverage.)
c. How do the average widths compare?
d. Is the normal interval's coverage close to 95%? If so, what does that tell you about whether the Markov interval's extra width is "worth it" here?
:::

::: {#ans-exr-markov-vs-normal-intervals .callout-answer .callout-note title="Template"}

*(your answer here)*

:::


# Chebyshev's Inequality

You've just derived a special case of Chebyshev's inequality. The general statement: for any random variable $X$ with finite variance,
$$
P(|X - \E[X]| \ge t) \le \frac{\Var[X]}{t^2}.
$$
When you applied Markov's inequality to $|\hat\theta - \theta|^2$ above, you got exactly this with $X = \hat\theta$.


# Biased Estimators

## Introduction

So far, we've exclusively talked about *unbiased estimators*. That is, estimators with the property that their expected value is equal to the estimation target.
$$
\hat\theta \qqtext{ is called unbiased if } \E\sb*{\hat{\theta}} = \theta.
$$
We say an estimator is *biased* if this isn't true. In this homework, we'll work with some biased estimators to get a sense of what bias does to our inference. The punchline: it messes up our interval estimates' coverage. We'll see exactly how much in next week's lecture.


## Using Prior Information

To get a sense of what bias means, let's consider a simple example of a biased estimator. Suppose we're estimating a population proportion $\theta$ using a sample $Y_1 \ldots Y_n$ drawn with replacement from a binary population. Instead of using the sample mean $\bar{Y}$, we use this estimator.
$$
\tilde{Y}_1 = \frac{1}{n+1}\cdot \cb*{\frac{1}{2} + \sum_{i=1}^{n}Y_i}
$$

What's going on here? We're mixing in a "prior observation" of $1/2$ with our sample. It's as if we'd seen one observation equal to $1/2$ before collecting any data, and we're averaging that with what we actually observe.


This is a specific example of a more general estimator that integrates information from a prior study. Maybe a real study or maybe a study we're imagining. Suppose we have $\nprior$ observations $\yprior_1 \ldots \yprior_{\nprior}$ from this study.^[Sometimes we call these *pseudo-observations*. This particular interpretation, in which we think of them as coming from a prior study, can be used to derive this estimator from Bayesian principles.]
$$
\thetaprior = \frac{1}{\nprior}\sum_{i=1}^{\nprior}\yprior_i.
$$

Averaging all our observations—from our current study and this prior one—gives us the following estimator.
$$
\begin{aligned}
\tilde{Y}_{\nprior}
&= \frac{1}{\nprior + n} \cb*{\sum_{i=1}^{\nprior}\yprior_i + \sum_{i=1}^{n} Y_i} \\
&= \frac{1}{\nprior + n}  \cb*{ \nprior\thetaprior + n\bar Y }
\end{aligned}
$$

We treat these prior observations, and therefore their mean $\thetaprior$, as deterministic.
The simple example we started with, $\tilde{Y}_1$, is a special case where we have a single prior observation $\yprior_1 = \frac{1}{2}$.

## Visualizing the Impact of Prior Information

Let's get a sense of how using prior information like this impacts our inference. The plot below shows the sampling distributions of three estimators of a population proportion $\theta$ at three sample sizes $n$: 10, 40, and 160. These are the estimators.

1. The estimator $\hat\theta_1 = \tilde{Y}_{\nprior}$ for $\thetaprior=3/4$ and $\nprior=10$ prior observations.
2. The estimator $\hat\theta_2 = \tilde{Y}_{n}$ for $\thetaprior=3/4$ and $\nprior=n$ prior observations. The bigger the sample, the more prior observations we use.
3. The sample mean, $\hat\theta_3 = \bar Y$.

As usual, the estimation target $\theta$ is indicated by a [green line]{.target}, the sampling distribution's mean by a solid [blue line]{.blue}, and its mean plus and minus two standard deviations by dotted [blue lines]{.blue}.

::: {#fig-3-sampling-distributions}
```{r}
#| label: plot-prior-information-grid
#| echo: false
#| fig-width: 10
#| fig-asp: 1
#| fig-alt: "A three-by-three grid of sampling distributions compares three estimators at sample sizes 10, 40, and 160. Green lines mark the target; blue solid and dotted lines mark each estimator's mean and two-standard-deviation range."
set.seed(1234)
n.prior = 10
theta.prior = .75
theta = .5

error = 2 * rnorm(100000)
bias = list(a = function(n) 0,
            b = function(n) n.prior * (theta.prior - theta) / (n.prior + n),
            c = function(n) n * (theta.prior-theta) / (n + n))                # n.prior = n
variance = list(a = function(n) theta*(1-theta)/n,
                b = function(n) theta*(1-theta) * n / (n.prior + n)^2,
                c = function(n) theta*(1-theta) * n / (n + n)^2)


plotdata = expand_grid(estimator = names(bias), n = c(10, 40, 160)) %>%
    pmap(function(estimator, n) {
        thebias = bias[[estimator]](n)
        thevar = variance[[estimator]](n)
        data.frame(estimator = estimator, n = n, muhat = thebias + error * sqrt(thevar))
    }) |> bind_rows()

summaries = plotdata |>
    group_by(n, estimator) |>
    summarize(mean = mean(muhat), sd = sd(muhat))

breaks = seq(-2/3, 2/3, by = 1/3)
signs = case_when(
    breaks < 0 ~ "-",
    breaks == 0 ~ "",
    breaks > 0 ~ "+"
)
offsets = case_when(
    breaks != 0 ~ sprintf("%d/3", 3*abs(breaks)),
    breaks == 0 ~ ""
)
labels.abstract = sprintf("θ %s %s", signs, offsets)

ggplot(plotdata) +
    geom_vline(aes(xintercept = theta), data = data.frame(theta=0), color='green',alpha=.4, linewidth=1.5) +
    geom_histogram(aes(x = muhat, y = after_stat(density)), bins = 100, alpha = .2) +
    geom_vline(aes(xintercept = mean), data = summaries, alpha=.4, color='blue') +
    geom_vline(aes(xintercept = mean + 2 * sd), data = summaries, linetype = "dotted", color='blue',alpha=.4) +
    geom_vline(aes(xintercept = mean - 2 * sd), data = summaries, linetype = "dotted", color='blue',alpha=.4) +
    scale_x_continuous(breaks = breaks, labels = labels.abstract, limits = c(-3/4,3/4)) +
    theme(axis.text.x = element_text(angle = -90, vjust = 0.5, hjust = 1)) +
    xlab("") +
    ylab("") +
    facet_grid(n ~ estimator, labeller = labeller(n = label_both, estimator = label_both))
```

Sampling distributions of three estimators at three sample sizes.
:::


# Bias, Variance, and Mean Squared Error

We call an estimator *consistent* if, when our sample size $n$ is large enough, it converges to the thing we intend to estimate. One way to talk about this is *convergence in mean square*: $\E[(\hat\mu - \mu)^2] \to 0$ as $n \to \infty$. The quantity $\E[(\hat\mu - \mu)^2]$ is the estimator's *mean squared error*. It decomposes as
$$
\underset{\text{mean squared error}}{\E\sb*{(\hat\mu - \mu)^2}} =  \underset{\text{bias}^2}{\p*{\E[\hat\mu] - \mu}^2} + \underset{\text{variance}}{\Var\sb*{\hat\mu}}.
$$
This is called the *bias/variance decomposition*. Since $\text{bias}^2$ and $\text{variance}$ are both positive, an estimator is consistent if and only if both its bias and its variance go to zero.

::: {.callout-exercise}
Using the bias and standard deviation formulas you calculated above, write a formula for the mean squared error of $\tilde{Y}_{\nprior}$. At what value of $\nprior$ is the mean squared error minimized? Does the answer depend on $\theta$?
:::

::: {#ans-q5 .callout-answer .callout-note title="Template"}

*(your answer here)*

:::


# Comparing Interval Calibration Methods

Now let's compare bootstrap, normal approximation, and Markov-calibrated intervals in practice. We'll use our GA turnout data.

::: {.callout-exercise}
For each of two point estimators --- (1) the sample mean $\hat\theta = \bar{Y}$ and (2) the prior-observations estimator $\tilde{Y}_{100}$ with $\thetaprior = 1/2$ and $\nprior=100$ ---

1. Plot the bootstrap sampling distribution and a 95% confidence interval calibrated using it.
2. Plot an estimate of the normal approximation to the estimator's sampling distribution and a second 95% confidence interval calibrated using that.
3. Plot a third 95% confidence interval, this time calibrated using Markov's inequality.
4. Taking advantage of your knowledge of the population, plot the estimator's actual sampling distribution and calculate the coverage probability of your three interval estimators.

Referring to your plots, comment on the behavior of your six^[$2\ \text{point estimators} \times 3\ \text{interval calibration methods}=6\ \text{interval estimators}$] interval estimators. If you had to do something like this again, which would you choose? Why? Would it make a difference if your sample size were larger, e.g. $n=2500$ instead of $n=625$?
:::

::: {#ans-q6 .callout-answer .callout-note title="Template"}

*(your answer here)*

:::
