import numpy as np
from numpy.random import default_rng
rng = default_rng(42)
from scipy.stats import norm
import matplotlib as mpl
from matplotlib import pyplot as plt
import plotly.graph_objects as go

mpl.rcParams['font.size'] = 18
plt.rcParams['text.usetex'] = False

1

Let $X \sim \text{Exponential}(\beta)$. Find $P(|X - \mu_X | \ge k \sigma_X)$ for $k \ge 1$. Compare this to the bound you get from Chebyshev's inequality.

Solution:

We have

\begin{align*} P(|X - \mu_X | \ge k \sigma_X) &= P(X \ge \beta(k+1)) + P(X \le \beta(1-k)) \\ &= P(X \ge \beta(k+1)) \tag{$P(X \le 0) = 0$} \\ &= 1 - \int_0^{\beta(k+1)} \frac{1}{\beta} e ^{-x / \beta} \, dx \\ &= e ^ {-(k+1)} \end{align*}

meanwhile, Chebyshev's inequality, with $t = k\sigma$, yields:

$$P(|X - \mu_X | \ge k \sigma_X) \le \frac{1}{k^2}$$

nn = np.arange(1, 10, 1)
plt.plot(nn, np.exp(-(nn + 1)), 'o-', label='$e^{-(k+1)}$')
plt.plot(nn, 1 / np.power(nn, 2), 'o-', label='$1 / k^2$')
plt.legend()
plt.grid()
plt.show()

2

Let $X \sim \text{Poisson}(\lambda)$. Use Chebyshev's inequality to show that $P(X \ge 2 \lambda) \le \frac{1}{\lambda}$.

Solution:

We have

\begin{align*} P(X \ge 2 \lambda) &= P(X - \lambda \ge \lambda) \\ &= P(X - \lambda \ge \lambda) + P(\lambda - X \ge \lambda) \tag{$P(X \le 0) = 0$} \\ &= P(|X - \lambda| \ge \lambda) \\ &\le \frac{\lambda}{\lambda^2} \tag{Chebyshev's Inequality with $t = \lambda$} \\ &= \frac{1}{\lambda} \end{align*}

3

Let $X_1, \dots, X_n \sim \text{Bernoulli}(p)$ and $\bar{X}_n = n^{-1} \sum_{i=1}^n X_i$. Bound $P(|\bar{X}_n - p| > \epsilon)$ using Chebyshev's inequality and using Hoeffding's inequality. Show that, when $n$ is large, the bound from Hoeffding's inequality is smaller than the bound from Chebyshev's inequality.

Solution:

We have $\mu = E(\bar{X}_n) = p$ and $\sigma^2 = V(\bar{X}_n) = \frac{p(1-p)}{n}$.

The Chebyshev bound is

$$P(|\bar{X}_n - p| > \epsilon) \le \frac{p(1-p)}{n \epsilon ^ 2}$$

and the Hoeffding bound is

$$P(|\bar{X}_n - p| > \epsilon) \le 2 e ^ {-2 n \epsilon ^ 2}$$

Let $n \ge \frac{1}{p(1-p)\epsilon^2}$. Let $A = p(1-p)$ and $B = 2\epsilon^2$.

We have

\begin{align*} \frac{2}{A B} &\le n \\ \Rightarrow n &\le \frac{AB^2n^2}{2B} \\ &= \frac{A}{B} \left[ \frac{B^2n^2}{2}\right] \\ &\le \frac{A}{B} e^{Bn} \tag{Taylor Expansion of $e^{Bn}$} \\ \Rightarrow 2e^{-Bn} &\le \frac{A}{(B/2)n} \Rightarrow 2e^{-2n\epsilon^2} \le \frac{p(1-p)}{n \epsilon^2} \end{align*}

4

Let $X_1, \dots, X_n \sim \text{Bernoulli}(p)$.

(a) Let $\alpha > 0$ be fixed and define

$$\epsilon_n = \sqrt{ \frac{1}{2n} \log \left( \frac{2}{\alpha} \right) }$$

Let $\hat{p}_n = n^{-1} \sum_{i=1}^n X_i$. Define $C_n = (\hat{p}_n - \epsilon_n, \hat{p}_n + \epsilon_n)$. Use Hoeffding's inequality to show that

$$P(\text{$C_n$ contains $p$}) \ge 1 - \alpha.$$

In practice, we truncate the interval so it does not go below 0 or above 1.

Solution:

Observe:

\begin{align*} P(\text{$C_n$ contains $p$}) &= P(|\hat{p}_n - p| \le \epsilon_n) \\ &= 1 - P(|\hat{p}_n - p| \ge \epsilon_n) \\ &\ge 1 - 2\exp \left[ -2 n \epsilon_n^2 \right] \tag{Hoeffding's inequality} \\ &= 1 - \alpha \end{align*}

where the last inequality holds because \begin{align*} 2 \exp \left[ -2 n \epsilon_n^2 \right] &= 2 \exp \left[ -2 n \left(\sqrt{ \frac{1}{2n} \log \left( \frac{2}{\alpha} \right) } \right)^2 \right] \\ &= \alpha. \end{align*}

(b) (Computer Experiment.) Let's examine the properties of this confidence interval. Let $\alpha = 0.05$ and $p = 0.4$. Conduct a simulation study to see how often the interval contains $p$ (called the coverage). Do this for various values of $n$ between 1 and 10000. Plot the coverage versus $n$.

alpha = 0.05
p = 0.4

N = 10000
k = 5000 # num simulations per sample size

nn = np.arange(1, N + 1)
coverage = np.empty(nn.size)
for n in nn:
    eps = np.sqrt(1 / (2 * n) * np.log(2 / alpha))
    p_hats = rng.binomial(n, p, size=k) / n
    coverage[n-1] = np.sum(np.where(np.abs(p_hats - p) <= eps, 1, 0)) / k

plt.figure(figsize=(12, 8))
plt.plot(nn,coverage)
plt.xlabel('n')
plt.ylabel('Coverage')
plt.ylim(0.95, 1)
plt.grid()
plt.show()

The coverage is higher (99.5%) than what might be suggested by the bound (95%).

(c) Plot the length of the interval versus $n$. Suppose we want the length of the interval to be no more than $.05$. How large should $n$ be?

nn = np.arange(1, 4000)
epsilon = np.sqrt( (1 / (2 * nn)) * np.log(2 / alpha))
width = 2 * epsilon
plt.plot(nn, width, label='width')
plt.plot(nn, [0.05 for n in nn], label='threshold')
plt.xlabel('$n$')
plt.yscale('log')
plt.grid()
plt.legend()
plt.show()

The minimum value number of observations required for a width of the $1-\alpha$ confidence interval under a value $W$, $n^*$, is given by

$$n^* = \left\lceil \frac{2}{W^2} \log \left( \frac{2}{\alpha} \right) \right\rceil. $$

If $W=0.05$ and $\alpha=0.05$, then $n^* = 2952$.

5

Prove Mill's inequality, Theorem 4.7:

Let $Z \sim N(0,1)$. Then,

$$P(|Z| > t) \le \sqrt{\frac{2}{\pi}} \frac{e^{-t^2 / 2}}{t}.$$

Hint: Note that $P(|Z| > t) = 2P(Z > t)$. Now write out what $P(Z > t)$ means and note that $x / t > 1$ whenever $x > t$.

Solution:

\begin{align*} P(|Z| > t) &= 2P(Z > t) \tag{symmetry}\\ &= 2 \frac{1}{\sqrt{2 \pi}} \int_t^{\infty} e^{\frac{-x^2}{2}} \, dx \tag{definition}\\ &\le \sqrt{\frac{2}{\pi}} \int_t^\infty \frac{x}{t} e^{\frac{-x^2}{2}} \, dx \tag{$\frac{x}{t} \ge 1$ inside integral}\\ &= \sqrt{\frac{2}{\pi}} \frac{e^{-t^2 / 2}}{t} \tag{u-substitution} \end{align*}

6

Let $Z \sim N(0,1)$. Find $P(|Z| > t)$ and plot this as a function of $t$. From Markov's inequality, we have the bound $P(|Z| > t) \le \frac{E\left[|Z|^k\right]}{t^k}$ for any $k > 0$. Plot these bounds for $k = 1,2,3,4,5$ and compare them to the true value of $P(|Z| > t)$. Also, plot the bound from Mill's inequality.

Solution:

n = 10000
tt = np.arange(0.1, 5, 0.1)
kk = np.arange(1, 6, 1)
Z = rng.normal(loc=0, scale=1, size=n)
probs = np.empty(tt.size)
for i,t in enumerate(tt):
    probs[i] = np.sum(np.abs(Z) > t) / n
markov_bounds = np.empty((kk.size, tt.size))
for k in kk:
    Z_power = np.abs(Z) ** k
    for i, t in enumerate(tt):
        markov_bounds[k-1, i] = np.mean(Z_power) / t ** k
mill_bound = np.sqrt(2 / np.pi) * np.exp(-np.power(tt, 2) / 2) / tt

plt.figure(figsize=(8, 6))
plt.plot(tt, probs, label='True Probability')
for k in kk:
    plt.plot(tt, markov_bounds[k-1,:], label=f'Markov Bound (k={k})')
plt.plot(tt, mill_bound, label='Mill Bound')
plt.ylim(0,1)
plt.legend()
plt.grid()
plt.show()

7

Let $X_1, \dots, X_n \sim N(0,1)$. Bound $P(|\bar{X}_n| > t)$ using Mill's inequality, where $\bar{X}_n = n ^{-1} \sum_{i=1}^n X_i$. Compare to the Chebyshev bound.

Solution:

Observe that $\bar{X}_n \sim N(0, \frac{1}{n})$, and that $\sqrt{n}\bar{X}_n \sim N(0, 1)$. Thus, we have the Mill bound:

\begin{align*} P(|\bar{X}_n| > t) &= P\left(\left| Z\right| > \sqrt{n}t\right) \\ &= \sqrt{\frac{2}{\pi n}} \frac{e^{-nt^2 / 2}}{t} \tag{Mill's inequality} \end{align*}

and the (weaker) Chebyshev bound:

\begin{align*} P(|\bar{X}_n| > t) &\le 1 / (nt ^ 2) \end{align*}

t = 0.3
nn = np.arange(1, 21, 1)
 
chebyshev_bound = 1 / (t * nn)
mill_bound = np.sqrt(2 / (np.pi * nn)) * np.exp(- nn * t ** 2 / 2) / t

fig = go.Figure()
fig.add_trace(go.Scatter(x=nn, y=chebyshev_bound, name='Chebyshev Bound'))
fig.add_trace(go.Scatter(x=nn, y=mill_bound, name='Mill Bound'))
fig.update_layout(
    title=f"Standard Normal Sample Mean Bounds (t={t})",
    xaxis_title=r"n",
    yaxis_title=r"$P(|\bar{X}_n| > t)$",
    font=dict(
        family="Courier New, monospace",
        size=12,
        color="RebeccaPurple"
    ))
fig.show()