Week 7: Estimation, Part 2

Where we’ve been

Where we’ve been

• Under i.i.d random sampling from the target population, we’ve established that…

Where we’ve been

This week

• Convergence in distribution
  • As sample size \(n\) becomes large, is the sampling distribution well-approximated by a known distribution?
  • When a central limit theorem applies, the sampling distribution converges in distribution to a normal!
• Plug-in estimation
  • What if we want to learn other things about the population distribution?
  • A population variance? A population covariance?
  • Like the mean, these are all functionals of the population CDF
  • Under i.i.d. sampling, plugging in the sample CDF usually gives a good estimator!
• Density estimation
  • Histograms and binning
  • Kernel estimators (weighted averages w/ weights determined by closeness)

Big-\(O\) and little-\(o\) notation

• Sometimes we’ll want to use shorthand to characterize how fast a term shrinks (or grows)

• Little-\(o\) notation: \(a_n = o(b_n)\) means \(\dfrac{a_n}{b_n} \to 0\) as \(n \to \infty\)

  • \(a_n\) vanishes faster than \(b_n\), e.g. \(\frac{1}{n^2} = o\left(\frac{1}{n}\right)\)

• Big-\(O\) notation: \(a_n = O(b_n)\) means \(\left|\dfrac{a_n}{b_n}\right|\) stays bounded as \(n \to \infty\)

  • \(a_n\) vanishes no slower than \(b_n\), e.g. \(\Var(\bar{Y}_n) = \frac{\sigma^2}{n} = O\left(\frac{1}{n}\right)\)

• So \(o(1)\) is anything going to \(0\), and \(O(1)\) is anything that stays bounded

Big-\(O_p\) and little-\(o_p\) notation

• For random variables, the same shorthand is stated in terms of probabilities

• Little-\(o_p\) notation denotes convergence (in probability) to zero

  • \(Y_n = o_p(1)\) means \(\Pr(|Y_n| \ge \epsilon) \to 0\) as \(n \to \infty\), for every \(\epsilon > 0\)
  • e.g. by the WLLN, we can write \(\bar{Y}_n = \mu + o_p(1)\)

• Big-\(O_p\) notation denotes boundedness (in probability) of a sequence as \(n \to \infty\)

  • \(Y_n = O_p(1)\) means that for every \(\epsilon > 0\) there’s an \(M\) with \(\Pr(|Y_n| > M) < \epsilon\) for all sufficiently large \(n\)

• The rest of the orders are defined by rescaling

  • So \(Y_n = O_p(a_n)\) means that \(\frac{Y_n}{a_n} = O_p(1)\), and likewise for \(o_p\)

Big-\(O_p\) and little-\(o_p\) notation

• Useful properties
  • \(o_p(1) + o_p(1) = o_p(1)\)
  • \(o_p(1) + O_p(1) = O_p(1)\)
  • \(O_p(1)o_p(1) = o_p(1)\)
  • If a sequence is \(o_p(1)\), it is \(O_p(1)\) but not vice-versa
• Most common use: convergence rates.
  • How quickly does \(\bar{Y}_n\) converges to \(\mu\)?
  • It turns out \(\bar{Y}_n - \mu = O_p(1/\sqrt{n})\): the error shrinks at rate \(1/\sqrt{n}\)

Convergence in distribution

• By the Weak Law of Large Numbers (WLLN) \(\bar{Y}_n \xrightarrow{p} \mu\)
  • Tells us how the sampling distribution collapses to the point we want…
  • …but says nothing about the shape of the sampling distribution for any given \(n\)
• To approximate the sampling distribution in a large (but finite) sample, we need a notion of convergence to a random variable
  • Does the distribution of \(bar{Y}_n\) get close to some known distribution as \(n\) gets large?
• Convergence in distribution
  • A sequence of random variables \(X_n\) with CDFs \(F_n\) converges in distribution to a random variable \(X\) with CDF \(F\) if \[\lim_{n \to \infty} F_n(x) = F(x)\] at every \(x\) where \(F\) is continuous
• We’ll denote this as \(X_n \xrightarrow{d} X\).

Convergence in distribution

• We want to approximate the sampling distribution of \(\bar{Y}_n\). Does \(\bar{Y}_n\) itself converge in distribution?

  • Convergence in probability implies convergence in distribution
  • So we end up with a useless result: \(\bar{Y}_n \xrightarrow{p} \mu\), so \(\bar{Y}_n \xrightarrow{d} \mu\)
  • And recentering still gives \(\bar{Y}_n - \mu \xrightarrow{d} 0\)

• We need to analyze a sequence that does not collapse to a single point

• We’ll rescale \(\bar{Y}_n - \mu\) by some factor that guarantees the variance converges to a non-zero constant \[\Var\big(n^r(\bar{Y}_n - \mu)\big) = n^{2r} \cdot \frac{\sigma^2}{n} = \sigma^2 n^{2r - 1} \to \begin{cases} 0 & r < 1/2 \\ \sigma^2 & r = 1/2 \\ \infty & r > 1/2 \end{cases}\]

• So scaling the error by \(\sqrt{n}\) keeps the variance stable:

  • This is what we’ll analyze: \(\sqrt{n}(\bar{Y}_n - \mu)\)
  • Its mean (\(0\)) and variance (\(\sigma^2\)) never change with \(n\). What about its shape?

The Central Limit Theorem

The Central Limit Theorem

• Central Limit Theorem: Let \(Y_1, Y_2, \dotsc, Y_n\) be i.i.d. with \(\mathbb{E}[Y_i] = \mu\) and \(\Var(Y_i) = \sigma^2 < \infty\). Then \[\sqrt{n}(\bar{Y}_n - \mu) \xrightarrow{d} \mathcal{N}(0, \sigma^2)\]

• Equivalently, the standardized sample mean \(\dfrac{\bar{Y}_n - \mu}{\sigma/\sqrt{n}} \xrightarrow{d} \mathcal{N}(0, 1)\)

• In the limit, the only things that matter about \(Y_i\) after averaging are its mean and variance

  • Whatever the shape of \(Y_i\), the averages become (roughly) normal.

• We’ll use this to say that for large \(n\), the sample mean is approximately normal \[\bar{Y}_n \overset{a}{\sim} \mathcal{N}\left(\mu, \frac{\sigma^2}{n}\right)\]

• Let’s prove it!

Moment generating functions

• Our proof requires a new concept: a moment generating function

  • This is a useful mathematical trick which lets us analyze an entire distribution by studying just one function.
  • Then we can apply the tools of calculus - specifically Taylor approximations and limits of functions

• The moment generating function (MGF) of a random variable \(X\) is \[M_X(t) = \mathbb{E}\left[e^{tX}\right]\] so long as it’s finite for all \(t\) in some interval around \(0\)

• Taking the derivative and evaluating at \(t = 0\) generates the moments:

  • \(M_X(0) = 1\)
  • \(M_X'(0) = \mathbb{E}[X]\)
  • \(M_X''(0) = \mathbb{E}[X^2]\), …

Moment generating functions

• Properties that we need:
  1. Uniqueness: Two random variables with the same MGF have the same distribution
  2. Continuity: If \(M_{X_n}(t) \to M_X(t)\) for every \(t\) near \(0\), then \(X_n \xrightarrow{d} X\)
  3. Independent sums: If \(X_1 \indep X_2\), then \(M_{X_1 + X_2}(t) = \mathbb{E}\left[e^{tX_1}e^{tX_2}\right] = M_{X_1}(t)M_{X_2}(t)\)
  4. The MGF of the Normal: If \(X \sim \mathcal{N}(0, \sigma^2)\), then \(M_X(t) = e^{\sigma^2 t^2/2}\)
• Our goal: Show that the MGF of \(\sqrt{n}(\bar{Y}_n - \mu)\) converges to \(e^{\sigma^2 t^2/2}\)

Proving the Central Limit Theorem

• Define \(W_i = Y_i - \mu\), the centered draw, and let \(M(t) = \mathbb{E}\left[e^{tW_i}\right]\) be its MGF. Then \[M(0) = 1, \qquad M'(0) = \mathbb{E}[W_i] = 0, \qquad M''(0) = \mathbb{E}[W_i^2] = \sigma^2\]

• Writing the scaled sum of the \(W_i\) \[\sqrt{n}(\bar{Y}_n - \mu) = \frac{1}{\sqrt{n}}\sum_{i=1}^n W_i\]

• Call the MGF of this scaled error \(M_n(t)\). We first want to write \(M_n(t)\) in terms of \(M(t)\)

• We’ll assume the MGF \(M\) exists

  • The general proof uses the characteristic function \(\mathbb{E}[e^{itW_i}]\), which always exists.

Proving the Central Limit Theorem

• First, turn the sum in the exponent into a product \[\begin{aligned} M_n(t) &= \mathbb{E}\left[e^{t\sqrt{n}(\bar{Y}_n - \mu)}\right] && \small\color{#777777}{\text{definition of the MGF}} \\ &= \mathbb{E}\left[\exp\left(\tfrac{t}{\sqrt{n}}\textstyle\sum_{i=1}^n W_i\right)\right] && \small\color{#777777}{\text{scaled error = scaled sum}} \\ &= \mathbb{E}\left[\exp\left(\textstyle\sum_{i=1}^n \tfrac{t}{\sqrt{n}} W_i\right)\right] && \small\color{#777777}{\text{constant inside the sum}} \\ &= \mathbb{E}\left[\textstyle\prod_{i=1}^n \exp\left(\tfrac{t}{\sqrt{n}} W_i\right)\right] && \small\color{#777777}{e^{a + b} = e^a e^b} \end{aligned}\]

• Then use the i.i.d. properties \[\begin{aligned} \mathbb{E}\left[\textstyle\prod_{i=1}^n \exp\left(\tfrac{t}{\sqrt{n}} W_i\right)\right] &= \textstyle\prod_{i=1}^n \mathbb{E}\left[\exp\left(\tfrac{t}{\sqrt{n}} W_i\right)\right] && \small\color{#777777}{\text{independent}} \\ &= \textstyle\prod_{i=1}^n M\left(\tfrac{t}{\sqrt{n}}\right) && \small\color{#777777}{\text{identical: each $W_i$ has MGF $M$}} \\ &= \left[M\left(\tfrac{t}{\sqrt{n}}\right)\right]^n && \small\color{#777777}{\text{$n$ identical factors}} \end{aligned}\]

Proving the Central Limit Theorem

• Take a second-order Taylor expansion of \(M(t)\) around \(t = 0\) \[M(t) = M(0) + t\,M'(0) + \frac{t^2}{2}M''(0) + o(t^2) = 1 + 0 + \frac{\sigma^2 t^2}{2} + o(t^2)\] \(o(t^2)\) is the remainder - we use little-o notation to show that it goes to \(0\) faster than \(t^2\)

• Evaluate it at \(t/\sqrt{n}\), which goes to \(0\) as \(n\) grows \[M_n(t) = \left[M\left(\frac{t}{\sqrt{n}}\right)\right]^n = \left[1 + \frac{\sigma^2 t^2}{2n} + o\left(\frac{1}{n}\right)\right]^n\]

• Factor so this looks like the familiar compound interest formula. \[M_n(t) = \left(1 + \frac{a_n}{n}\right)^n, \qquad a_n = \frac{\sigma^2 t^2}{2} + n \cdot o\left(\frac{1}{n}\right)\]

• The remainder is \(o\left(\frac{1}{n}\right)\) and shrinks faster than \(\frac{1}{n}\), so even \(n\) times it goes to \(0\).

  • So \(a_n \to a = \frac{\sigma^2 t^2}{2}\)

Proving the Central Limit Theorem

• Since \(a_n \to a\), for large \(n\) our \(M_n(t)\) behaves just like compound interest at the fixed rate \(a\), which is one our definitions of \(e\). \[M_n(t) = \left(1 + \frac{a_n}{n}\right)^n \approx \left(1 + \frac{a}{n}\right)^n \to e^{a} = e^{\sigma^2 t^2/2}\]

• That’s the MGF of a \(\mathcal{N}(0, \sigma^2)\), so \(\sqrt{n}(\bar{Y}_n - \mu) \xrightarrow{d} \mathcal{N}(0, \sigma^2)\)

Many central limit theorems

• What we proved is the Lindeberg–Lévy CLT. There are many other CLTs, each under different data-generating processes

• The theme: a sum of many small, not-too-dependent pieces is approximately normal.
  • What changes is the sampling variance

Properties of convergence in distribution

• Consider \(X_n \xrightarrow{d} X\), \(Y_n \xrightarrow{p} c\) for a constant \(c\), and a continuous function \(g()\)

• By the continuous mapping theorem, \(g(X_n) \xrightarrow{d} g(X)\)

• By Slutsky’s theorem

1. \(X_n + Y_n \xrightarrow{d} X + c\)
2. \(X_nY_n \xrightarrow{d} cX\)
3. \(X_n/Y_n \xrightarrow{d} X/c\) if \(c \neq 0\)

• Careful: \(X_n \xrightarrow{d} X\) and \(Y_n \xrightarrow{d} Y\) does not imply \(X_n + Y_n \xrightarrow{d} X + Y\) in general

• Why this matters: we don’t know \(\sigma\), but if \(\hat{\sigma} \xrightarrow{p} \sigma\) (more on this later), we can plug it in and still know the asymptotic distribution. \[\frac{\bar{Y}_n - \mu}{\hat{\sigma}/\sqrt{n}} = \underbrace{\frac{1}{\hat{\sigma}}}_{\xrightarrow{p} 1/\sigma} \cdot \underbrace{\sqrt{n}(\bar{Y}_n - \mu)}_{\xrightarrow{d} \mathcal{N}(0, \sigma^2)} \xrightarrow{d} \frac{1}{\sigma}\mathcal{N}(0, \sigma^2) = \mathcal{N}(0, 1)\]

Averaging gets rid of everything but mean and variance

What we’ve accomplished

• Under i.i.d. sampling with a finite variance, the CLT fills in the question mark: in large samples, the sampling distribution of the sample mean is approximately normal
• How large is “large”? It depends on the distribution of \(Y_i\).

Plug-in estimators

• So far, we’ve estimated one feature of the population: the mean \(\mu = \mathbb{E}[Y_i]\), with the sample mean \(\bar{Y}_n\)

• But we often want to learn other features of the population distribution

  • The variance \(\sigma^2\) (which we need for standard errors)
  • The median, or some other quantile
  • The share of the population above some threshold, \(\Pr(Y_i > c)\)

• Informally: Write down the feature of the population we care about, then compute its sample analog

  • Expected value \(\rightarrow\) sample mean.
  • Variance \(\rightarrow\) sample variance.

• Each of these features is a function of the CDF \(F\)

  • Estimate \(F\) with its sample analog, then plug it in

The empirical CDF

• By the fundamental bridge, the CDF is an expectation of an indicator random variable

  \[F(y) = \Pr(Y_i \le y) = \mathbb{E}\left[I(Y_i \le y)\right]\]

• Its sample analog is the empirical CDF: replace the expectation with a sample mean \[\hat{F}_n(y) = \frac{1}{n}\sum_{i=1}^n I(Y_i \le y) = \text{share of observations at or below } y\]

• A step function that jumps by \(\frac{1}{n}\) at each observation, even when \(Y_i\) is continuous

The empirical CDF

• For any given value of \(y\), the indicators \(Z_i = I(Y_i \le y)\) are i.i.d. Bernoulli with mean \(\Pr(Y_i \le y) = F(y)\)

• So \(\hat{F}_n(y) = \bar{Z}_n\) is just a sample mean, and everything we know about sample means applies

  • Unbiased: \(\mathbb{E}\left[\hat{F}_n(y)\right] = F(y)\)
  • Sampling variance: \(\Var\left(\hat{F}_n(y)\right) = \dfrac{F(y)\big(1 - F(y)\big)}{n}\)
  • Consistent (WLLN): \(\hat{F}_n(y) \xrightarrow{p} F(y)\)
  • Asymptotically normal (CLT): \(\sqrt{n}\left(\hat{F}_n(y) - F(y)\right) \xrightarrow{d} \mathcal{N}\big(0, F(y)(1 - F(y))\big)\)

• Even stronger, the Glivenko–Cantelli theorem shows that the largest gap over all \(y\) also goes to zero (with probability 1) \[\sup_y \left|\hat{F}_n(y) - F(y)\right| \to 0\]

  • The whole curve converges

The plug-in principle

• A statistical functional \(T\) maps a CDF to a number. Many estimands can be written as \(\theta = T(F)\)

• Notation: \(\int g(y)\,dF(y)\) means \(\sum_y g(y) f(y)\) if \(Y_i\) is discrete, or \(\int g(y) f(y)\,dy\) if it’s continuous.

  • It’s \(\mathbb{E}[g(Y_i)]\)

• Examples

  • Mean: \(\mu = T_{\mu}(F) = \int y\,dF(y)\)
  • Variance: \(\sigma^2 = T_{\sigma^2}(F) = \int y^2\,dF(y) - \left(\int y\,dF(y)\right)^2\)
  • Median: \(T_{\text{med}}(F) = F^{-1}(1/2)\)
  • Tail probability: \(\Pr(Y_i > c) = 1 - F(c)\)

• The plug-in estimator of \(\theta = T(F)\) applies the same functional to the empirical CDF \[\hat{\theta} = T(\hat{F}_n)\]

The plug-in principle

• \(\hat{F}_n\) is the CDF of a discrete distribution that puts mass \(\frac{1}{n}\) on each observation, so for any \(g\) \[\int g(y)\,d\hat{F}_n(y) = \sum_{i=1}^n g(Y_i) \cdot \frac{1}{n} = \frac{1}{n}\sum_{i=1}^n g(Y_i)\]

  • Every expectation becomes a sample mean

• Mean: \(T_{\mu}(\hat{F}_n) = \frac{1}{n}\sum_{i=1}^n Y_i = \bar{Y}_n\).

• Variance: \(T_{\sigma^2}(\hat{F}_n) = \frac{1}{n}\sum_{i=1}^n Y_i^2 - \bar{Y}_n^2 = \frac{1}{n}\sum_{i=1}^n (Y_i - \bar{Y}_n)^2\), the plug-in sample variance

• Median: \(\hat{F}_n^{-1}(1/2)\), the sample median.

• Tail probability: \(1 - \hat{F}_n(c)\), the share of observations above \(c\)

Bias of the plug-in variance

• The population variance \[\sigma^2 = \mathbb{E}[Y_i^2] - \mu^2 \]

• We can re-write our plug-in variance estimator the same way.

  \[\begin{aligned} \hat{\sigma}^2_n &= \frac{1}{n}\sum_{i=1}^n (Y_i - \bar{Y}_n)^2 \\ &= \frac{1}{n}\sum_{i=1}^n \left(Y_i^2 - 2Y_i\bar{Y}_n + \bar{Y}_n^2\right) && \small\color{#777777}{\text{expand the square}} \\ &= \frac{1}{n}\sum_{i=1}^n Y_i^2 - 2\bar{Y}_n \cdot \frac{1}{n}\sum_{i=1}^n Y_i + \frac{1}{n}\sum_{i=1}^n \bar{Y}_n^2 && \small\color{#777777}{\text{split the sum; } \bar{Y}_n \text{ doesn't depend on } i} \\ &= \frac{1}{n}\sum_{i=1}^n Y_i^2 - 2\bar{Y}_n \cdot \bar{Y}_n + \bar{Y}_n^2 && \small\color{#777777}{\tfrac{1}{n}\textstyle\sum_i Y_i = \bar{Y}_n \text{, and } n \text{ copies of } \bar{Y}_n^2} \\ &= \frac{1}{n}\sum_{i=1}^n Y_i^2 - \bar{Y}_n^2 \end{aligned}\]

Bias of the plug-in variance

• Is it unbiased? Take the expectation of each piece. First, the average of the squares \[\begin{aligned} \mathbb{E}\left[\frac{1}{n}\sum_{i=1}^n Y_i^2\right] &= \frac{1}{n}\sum_{i=1}^n \mathbb{E}\left[Y_i^2\right] && \small\color{#777777}{\text{linearity of expectations}} \\ &= \frac{1}{n} \cdot n\,\mathbb{E}\left[Y_i^2\right] = \mathbb{E}\left[Y_i^2\right] && \small\color{#777777}{\text{identically distributed}} \end{aligned}\]

• Then the square of the average. Rearrange \(\Var(X) = \mathbb{E}[X^2] - \mathbb{E}[X]^2\) with \(X = \bar{Y}_n\) \[\begin{aligned} \mathbb{E}\left[\bar{Y}_n^2\right] &= \Var(\bar{Y}_n) + \mathbb{E}[\bar{Y}_n]^2 \\ &= \frac{\sigma^2}{n} + \mu^2 && \small\color{#777777}{\text{last week: } \mathbb{E}[\bar{Y}_n] = \mu,\ \Var(\bar{Y}_n) = \sigma^2/n} \end{aligned}\]

  • \(\mathbb{E}[\bar{Y}_n^2] \neq \mathbb{E}[\bar{Y}_n]^2\): the sample mean’s own sampling variance appears!

Bias of the plug-in variance

• Combining the components \[\begin{aligned} \mathbb{E}\left[\hat{\sigma}^2_n\right] &= \mathbb{E}\left[\frac{1}{n}\sum_{i=1}^n Y_i^2\right] - \mathbb{E}\left[\bar{Y}_n^2\right] && \small\color{#777777}{\text{linearity of expectations}} \\ &= \mathbb{E}\left[Y_i^2\right] - \left(\frac{\sigma^2}{n} + \mu^2\right) && \small\color{#777777}{\text{the two pieces}} \\ &= \left(\mathbb{E}\left[Y_i^2\right] - \mu^2\right) - \frac{\sigma^2}{n} && \small\color{#777777}{\text{regroup}} \\ &= \sigma^2 - \frac{\sigma^2}{n} && \small\color{#777777}{\text{definition of } \sigma^2} \\ &= \frac{n - 1}{n}\sigma^2 \end{aligned}\]

• So it’s biased downward: on average, it’s too small by \(\frac{\sigma^2}{n}\)

  • Intuition: Think about computing the plug-in variance with just a single observation!

Fixing the bias

• How big is the bias? \(\mathbb{E}\left[\hat{\sigma}^2_n\right] = \frac{n-1}{n}\sigma^2\) is too small by a fraction \(\frac{1}{n}\)

• Rescale by \(\frac{n}{n-1}\) to get the unbiased sample variance estimator (R’s var()) \[S^2_n = \frac{n}{n-1}\hat{\sigma}^2_n = \frac{1}{n-1}\sum_{i=1}^n (Y_i - \bar{Y}_n)^2 \qquad \mathbb{E}\left[S^2_n\right] = \frac{n}{n-1} \cdot \frac{n-1}{n}\sigma^2 = \sigma^2\]

• Dividing by \(n - 1\) instead of \(n\) eliminates the bias (“Bessel’s correction”)

• But the bias vanishes as \(n\) grows. Does \(\hat{\sigma}^2_n\) still get the right answer in large samples?

Consistency of the plug-in variance

• The \(Y_i^2\) are i.i.d. with mean \(\mathbb{E}[Y_i^2]\), so by the WLLN \[\frac{1}{n}\sum_{i=1}^n Y_i^2 \xrightarrow{p} \mathbb{E}[Y_i^2]\]

• By the WLLN, \(\bar{Y}_n \xrightarrow{p} \mu\). A limit of a product is the product of the limits, so \(\bar{Y}_n^2 = \bar{Y}_n \cdot \bar{Y}_n \xrightarrow{p} \mu^2\)

• A sum (or difference) converges to the sum of the limits, so \[\hat{\sigma}^2_n = \frac{1}{n}\sum_{i=1}^n Y_i^2 - \bar{Y}_n^2 \xrightarrow{p} \mathbb{E}[Y_i^2] - \mu^2 = \sigma^2\]

• Biased, but consistent

Consistency of plug-in estimators

• In general: \(\hat{F}_n \to F\). So if small changes in \(F\) only cause small changes in \(T(F)\), then \(T(\hat{F}_n) \to T(F)\)

  • In other words, \(T\) needs to be continuous (in a particular sense)

• Under some additional regularity conditions, plug-in estimators are also asymptotically normal

  • This involves a kind of differentiability of \(T\) with respect to the distribution \(F\).

• Most estimands we’ll work with are well behaved: means, variances, covariances, quantiles (where \(F\) isn’t flat), regression coefficients

• When does it fail? When \(T\) is too sensitive to small changes in \(F\)

  • The density \(f(y) = F'(y)\): \(\hat{F}_n\) is a step function and fails continuity
  • Plugging \(\hat{F}_n\) straight doesn’t work. We’ll need to smooth it.

Running Example: 2024 GSS

• Estimate features of the age distribution of U.S. adults from the 3,859 GSS respondents, treating them as an i.i.d. sample

• With \(n\) this large, the \(\frac{n}{n-1}\) correction barely matters: \(S^2_n = 312.0\)
• And the standard error of the mean: \(\widehat{SE}(\bar{Y}_n) = S_n / \sqrt{n} = 0.284\)

Running Example in R

library(tidyverse)

gss <- read_csv("../week6/data/gss2024.csv") |>
  filter(!is.na(age))       # the GSS top-codes age at 89

Fhat <- ecdf(gss$age)       # the empirical CDF, as an R function

age_est <- gss |>
  summarize(n = n(),
            mean = mean(age),
            var_plugin = mean((age - mean(age))^2),
            var_unbiased = var(age),
            median = median(age),
            share65 = 1 - Fhat(64))
age_est

# A tibble: 1 × 6
      n  mean var_plugin var_unbiased median share65
  <int> <dbl>      <dbl>        <dbl>  <dbl>   <dbl>
1  3859  49.0       312.         312.     48   0.236

Density estimation

• The plug-in principle can’t give us the density \(f(y)\) directly

  • But the density is often the best way to visualize the shape of a distribution

• Recall what a density is: the share of the population in a small window around \(y\), divided by the window’s width \[f(y) \approx \frac{F(y + h/2) - F(y - h/2)}{h} \qquad \text{for small } h\]

• Plug in \(\hat{F}_n\): the share of the sample within \(h/2\) of \(y\), divided by \(h\) \[\hat{f}(y) = \frac{\hat{F}_n(y + h/2) - \hat{F}_n(y - h/2)}{h} = \frac{\#\{i : |Y_i - y| \le h/2\}}{nh}\]

• The width \(h\) is a tuning parameter

  • Too small and each window holds few observations (noisy)
  • Too large and we average over a region where \(f\) changes (biased)

Histograms

• A histogram splits the real number line into bins \(B_1, B_2, \dotsc\) of width \(h\). For \(y\) in bin \(B_j\) \[\hat{f}(y) = \frac{\#\{i : Y_i \in B_j\}}{nh}\]

  • Each bar’s area is the share of observations in its bin, so the areas sum to \(1\)

• Bias: the bar is flat but \(f\) might not be. Larger \(h\) \(\leadsto\) larger bias (depending on \(f^\prime(y)\))

• Variance: the count in a bin is binomial with \(n\) trials and success probability approximately \(hf(y)\) \[\Var\left(\hat{f}(y)\right) \approx \frac{n hf(y) (1 - hf(y))}{(nh)^2} = \frac{f(y)}{nh} (1 - hf(y)) \approx \frac{f(y)}{nh} \qquad \text{for small } h\]

• Consistent if the bins shrink, \(h \to 0\), but slowly enough that each still fills up, \(nh \to \infty\)

Choosing the bin width

• Pick \(h\) to minimize the integrated mean squared error across all \(y\). For small \(h\) (Scott 1979) \[\text{IMSE}(h) \approx \underbrace{\frac{h^2}{12}\int f'(y)^2\,dy}_{\text{bias}^2} + \underbrace{\frac{1}{nh}}_{\text{variance}} \quad \Longrightarrow \quad h^* = \left(\frac{6}{\int f'(y)^2\,dy}\right)^{1/3} n^{-1/3}\]
  • Bins shrink at a rate of \(n^{-1/3}\), so the number of bins grows like \(n^{1/3}\)
• \(h^*\) depends on the unknown \(f\) through \(\int f'(y)^2\,dy\). Heuristics plug in a reference distribution
  • Scott (normal reference distribution): \(h = 3.49\,\hat{\sigma}\,n^{-1/3}\) (Scott 1979)
  • Freedman–Diaconis: \(h = 2\,\text{IQR}\,n^{-1/3}\) (Freedman and Diaconis 1981)
  • Sturges: \(\lceil \log_2 n \rceil + 1\) bins (Sturges 1926). This is R’s default in hist() - tends to use too few bins in large samples

How many bins?

• Too narrow: mostly noise. Too wide: the shape gets smoothed away

Kernel density estimation

• A histogram’s estimate at \(y\) depends on where the bin edges happen to fall

• Instead, we could imagine centering a window on every point \(y\) (Rosenblatt 1956) \[\hat{f}(y) = \frac{1}{nh}\sum_{i=1}^n I\left(|Y_i - y| \le \tfrac{h}{2}\right)\]

• A weighted count: observations within \(h/2\) of \(y\) get weight \(1\). Everyone else gets nothing

  • But instead of weighting every observation in the window equally, we could let an observation’s weight fall off smoothly with its distance from \(y\)

Kernel density estimation

• The kernel density estimator (Rosenblatt 1956; Parzen 1962)

  \[\hat{f}_h(y) = \frac{1}{nh}\sum_{i=1}^n K_h(Y_i - y) \]

• Where the kernel is defined as:

  \[K_h(u) = K\left(\frac{u}{h}\right)\]

• The kernel \(K\) is a function that is symmetric, non-negative and integrates to \(1\):

  • e.g. uniform, triangular \(1 - |u|\) on \([-1, 1]\), Gaussian \(\phi(u)\)
  • The choice of kernel defines the shape of how the weights decay from \(y\)

• The bandwidth \(h\) defines how many observations get meaningful weight.

Kernel density estimation

• Each observation contributes a small bump \(\frac{1}{nh}K_h(\cdot - Y_i)\) (dashed). The estimate is their sum (solid)

Kernel density estimation

• Same observations with a triangular kernel instead of a gaussian/normal kernel

Running Example: 2024 GSS

Running Example in R

h_fd   <- 2 * IQR(gss$age) / nrow(gss)^(1/3)   # Freedman-Diaconis bin width
bw_rot <- bw.nrd0(gss$age)                     # Silverman's rule of thumb
c(bin_width = h_fd, bandwidth = bw_rot)

bin_width bandwidth 
 3.825249  3.047998 

fhat <- density(gss$age, bw = bw_rot)          # Gaussian kernel by default
approx(fhat$x, fhat$y, xout = c(30, 60))$y     # estimated density at ages 30 and 60

[1] 0.01749374 0.01663595

ggplot(gss, aes(age)) +
  geom_histogram(aes(y = after_stat(density)), binwidth = h_fd) +
  geom_density(bw = bw_rot)

Next week

• What do we do now that we have the sampling distribution of the sample mean?
• Interval estimators
  • Instead of a single point estimate can we give a range of values that characterizes our uncertainty?
  • Confidence interval: In some fraction (e.g. 95%) of samples, the bounds will cover the true parameter.
  • For a consistent and asymptotically normal estimator, this is just \(\hat{\mu} \pm 1.96 \times \widehat{SE(\hat{\mu})}\) (for a 95% CI).
• Hypothesis testing
  • A thought experiment: If the true population mean were \(\mu_0\), how likely would I be to observe a realization as extreme as the one in my sample?
  • This is a p-value
  • If the value is sufficiently unlikely…I may be willing to believe that the true population mean is not \(\mu_0\)