Симуляция доверительных интервалов методом Монте-Карло

import numpy as np
import scipy.stats as stats
import seaborn as sns
from matplotlib import pyplot as plt

Доверительные интервалы

Иногда не так важно узнать точное значение параметра — достаточно указать область, где он, вероятно, находится. Вернёмся к самому первому примеру с оценкой ежедневной прибыли будущего бизнеса продаж. Необязательно утверждать абсолютно точно, что «средняя дневная прибыль будет от \(1800\) до \(2200\)», но это значительно лучше, чем ошибочно сказать «средняя дневная прибыль равна \(4000\)».

Самый простой и естественный способ провести точечную оценку неизвестного математического ожидания \(E[X]\) — это, конечно, выборочное среднее \(\bar{X}\). Но использовать выборочное среднее как точную оценку \(\mu\) может быть рискованно из-за случайной природы выборки. Если повезёт — реализация \(\bar{X}\) окажется близкой к истинному значению, если нет — в выводах появится заметная ошибка.

Поэтому при использовании только точечной оценки мы, скорее всего, не попадаем точно в неизвестный параметр.

Но если будем использовать диапазон правдоподобных значений — доверительный интервал — то с высокой вероятностью «накрываем» неизвестный параметр.

За формулами и их подробными выводами см. лекции и другие материалы. В этом ноутбуке мы фокусируемся на важной не всегда очевидной интерпретации доверительных интервалов.

Так как \(L\) и \(U\), границы интервала, являются функциями от случайной выборки, они сами случайны.

Для каждой новой выборки реализации границ интервала будут разными.

Уровень доверия \((1-\alpha)\) означает, что среди всех возможных выборок примерно в \((1-\alpha)\cdot 100\%\) случаев интервал будет содержать неизвестный параметр.

1.a Доверительные интервалы для математического ожидания \(\mu\) при известной дисперсии

Продолжим пример с бизнесом. Пусть у нас есть выборка размера \(n\) из ежедневной прибыли аналогичного бизнеса на соседней улице. Также предположим, что прибыль имеет нормальное распределение \(X \sim \mathcal{N}(2000, 500^2)\).

unknown_var = stats.norm(2000, 500)
mu = unknown_var.mean()
var = unknown_var.var()
std = unknown_var.std()

n_draws = 100
sample_size = 10

alpha = 0.2
z_a = stats.norm(0, 1).ppf(1 - alpha / 2)  # (1-alpha) * 1000% доверительный интервал

result = np.empty(n_draws)
for i in range(n_draws):
    sample_realization = unknown_var.rvs(size=sample_size)
    result[i] = np.mean(sample_realization)

upper_bound = result + z_a * std / np.sqrt(sample_size)
lower_bound = result - z_a * std / np.sqrt(sample_size)

Сначала сделаем не очень много выборок, чтобы визуально увидеть: границы реализации доверительного интервала меняются при каждой новой выборке.

figure = plt.figure(figsize=(10, 5))
ax = figure.add_subplot()
ax.scatter(np.arange(n_draws), result, label=r'реализация выборочного среднего $\bar{x}$')
ax.errorbar(np.arange(n_draws), result, yerr=upper_bound - lower_bound, ls='none')
ax.axhline(mu, color='green', label = r'Истинное среднее, $\mu$')

plt.legend(loc='best')
figure.suptitle(r'Объем выборки = %d, истинные параметры распределения: $\mu = %s, \sigma = %s$'% (sample_size, mu, std))
plt.show()

Теперь соберем статистику. Сделаем много выборок и посчитаем, как часто доверительный интервал накрывает истинный параметр.

n_draws = 50000
n_correct_captures = 0
error = z_a * std / np.sqrt(sample_size)
for i in range(n_draws):
    sample_realization = unknown_var.rvs(size=sample_size)
    mean_point_estimate = np.mean(sample_realization)
    if (mu < mean_point_estimate + error ) and (mu > mean_point_estimate - error):
        n_correct_captures += 1

print(n_correct_captures / n_draws)

И действительно, доля покрытий близка к заявленному уровню доверия.

1.b Доверительные интервалы для истинной доли \(p\)

Минимальная теория (по мотивам лекции 18):

Пусть \(X_i \sim \text{Bernoulli}(p)\). Точечная оценка истинной доли — выборочная доля \[ \hat p = \frac{1}{n}\sum_{i=1}^n X_i. \] Для достаточно больших выборок (условно, \(n\hat p \gtrsim 5\), и \(n(1-\hat p) \gtrsim 5\), и \(n \geq 30\)), по ЦПТ используем нормальную аппроксимацию и получаем приближенный \((1-\alpha)\) доверительный интервал.

Теоретический результат: \[ p \in \left(\hat p - z_{\alpha/2}\sqrt{\frac{p(1-p)}{n}},\; \hat p + z_{\alpha/2}\sqrt{\frac{p(1-p)}{n}}\right). \]

На практике реализация \((1-\alpha)\) доверительного интервала:

\[ p \in \left(\tilde{p} - z_{\alpha/2}\sqrt{\frac{\tilde{p}(1-\tilde{p})}{n}},\; \tilde{p} + z_{\alpha/2}\sqrt{\frac{\tilde{p}(1-\tilde{p})}{n}}\right), \]

где \(\tilde{p}\) - реализация выборочной доли на реализации случайной выборки.

Как и в пункте 1.a, проверим покрытие интервала через Monte Carlo-симуляцию.

# 1.b: параметры симуляции
p_true = 0.39
bern = stats.bernoulli(p_true)

n_draws = 100
sample_size = 80

alpha = 0.15
z_a = stats.norm(0, 1).ppf(1 - alpha / 2)  # (1-alpha) * 1000% доверительный интервал

tilde_p = np.empty(n_draws)
lower_bound = np.empty(n_draws)
upper_bound = np.empty(n_draws)

for i in range(n_draws):
    sample_realization = bern.rvs(size=sample_size)
    tilde_p[i] = np.mean(sample_realization)
    error = z_a * np.sqrt(tilde_p[i] * (1 - tilde_p[i]) / sample_size)
    lower_bound[i] = tilde_p[i] - error
    upper_bound[i] = tilde_p[i] + error

Сначала сделаем немного выборок, чтобы визуально увидеть: границы интервала меняются для каждой новой выборки (теперь для доли \(p\)).

figure = plt.figure(figsize=(10, 5))
ax = figure.add_subplot()
ax.scatter(np.arange(n_draws), tilde_p, label=r'реализация выборочной доли $\tilde{p}$')
ax.errorbar(np.arange(n_draws), tilde_p, yerr=upper_bound - lower_bound, ls='none')
ax.axhline(p_true, color='green', label=r'истинная доля, $p$')

plt.legend(loc='best')
figure.suptitle(r'Объем выборки = %d, параметр Бернулли: $p = %.2f$' % (sample_size, p_true))
plt.show()

Теперь соберем статистику. Сделаем много выборок и посчитаем, как часто доверительный интервал накрывает истинную долю \(p\).

n_draws = 50000
n_correct_captures = 0

for i in range(n_draws):
    sample_realization = bern.rvs(size=sample_size)
    tilde_p = np.mean(sample_realization)
    error = z_a * np.sqrt(tilde_p * (1 - tilde_p) / sample_size)
    if (p_true < tilde_p + error) and (p_true > tilde_p - error):
        n_correct_captures += 1

print(n_correct_captures / n_draws)

При достаточно большом объеме выборки эмпирическое покрытие близко к номинальному уровню доверия.

1.c Доверительные интервалы для математического ожидания \(\mu \equiv E[X]\) при неизвестной дисперсии

Минимальная теория (по мотивам лекции 19):

Если дисперсия генеральной совокупности неизвестна, для нормальной выборки используем статистику Стьюдента \[ \frac{\bar X - \mu}{S/\sqrt{n}} \sim t_{n-1}, \quad S^2 = \frac{1}{n-1}\sum_{i=1}^n (X_i-\bar X)^2. \]

Тогда теоретический \((1-\alpha)\) доверительный интервал для математического ожидания: \[ \mu \in \left(\bar X - t_{n-1,\alpha/2}\frac{S}{\sqrt{n}},\; \bar X + t_{n-1,\alpha/2}\frac{S}{\sqrt{n}}\right). \]

На практике реализация такого \((1-\alpha)\) доверительного интервала: \[ \mu \in \left(\bar x - t_{n-1,\alpha/2}\frac{s}{\sqrt{n}},\; \bar x + t_{n-1,\alpha/2}\frac{s}{\sqrt{n}}\right), \]

где \(\bar{x}\) и \(s\) - реализации выборочного среднего и выборочного стандартного отклонения соответственно.

Проверим покрытие интервала Monte Carlo-симуляцией по логике пункта 1.a.

# 1.c: параметры симуляции
unknown_var = stats.norm(2000, 500)
mu_true = unknown_var.mean()

n_draws = 100
sample_size = 10
alpha = 0.15
t_a = stats.t(df=sample_size - 1).ppf(1 - alpha / 2)

sample_means = np.empty(n_draws)
lower_bound = np.empty(n_draws)
upper_bound = np.empty(n_draws)

for i in range(n_draws):
    sample_realization = unknown_var.rvs(size=sample_size)
    mean_estimate = np.mean(sample_realization)
    s_estimate = np.std(sample_realization, ddof=1)
    error = t_a * s_estimate / np.sqrt(sample_size)

    sample_means[i] = mean_estimate
    lower_bound[i] = mean_estimate - error
    upper_bound[i] = mean_estimate + error

Сначала сделаем немного выборок, чтобы визуально увидеть: границы доверительного интервала меняются при каждой новой выборке (теперь при неизвестной дисперсии и \(t\)-критическом значении).

figure = plt.figure(figsize=(10, 5))
ax = figure.add_subplot()
ax.scatter(np.arange(n_draws), sample_means, label=r'реализация выборочного среднего $\bar{x}$')
ax.errorbar(np.arange(n_draws), sample_means, yerr=upper_bound - lower_bound, ls='none')
ax.axhline(mu_true, color='green', label=r'истинное среднее, $\mu$')

plt.legend(loc='best')
figure.suptitle(r'Объем выборки = %d, случай неизвестной дисперсии ($t$-интервал)' % sample_size)
plt.show()

Теперь соберем статистику. Сделаем много выборок и посчитаем, как часто \(t\)-доверительный интервал накрывает истинное среднее \(\mu\).

n_draws = 50000
n_correct_captures = 0

for i in range(n_draws):
    sample_realization = unknown_var.rvs(size=sample_size)
    mean_point_estimate = np.mean(sample_realization)
    s_point_estimate = np.std(sample_realization, ddof=1)
    error = t_a * s_point_estimate / np.sqrt(sample_size)
    if (mu_true < mean_point_estimate + error) and (mu_true > mean_point_estimate - error):
        n_correct_captures += 1

print(n_correct_captures / n_draws)

Здесь эмпирическое покрытие тоже близко к номинальному уровню доверия.