<jupyter_start><jupyter_text>7. Доверительные интервалы. Пусть выборка $X_1, ..., X_n$ порождается некоторым распределением $F$. Пусть у распределения $F$ есть некий параметр $\theta$. Например это может быть вероятность успеха в схеме Бернулли, коэффициент эксцесса распределения $F$, или какой-нибудь другой функционал $T(F)$. Зададимся некоторым числом $\alpha$. Наша цель - написать правило $R(X_1, ..., X_n)$. Это правило будет по выборке выдавать 2 числа - начало и конец промежутка, в котором лежит $\theta$. Этот промежуток называют доверительным интервалом. Поскольку $R$ зависит от случайных величин, $R$ сама также является случайной величиной. Поэтому возникает ещё одно требование к правилу: для заданного в условии уровня $\alpha$ правило должно выдавать такой промежуток, что $\theta$ лежит в нём с вероятностью $1 - \alpha$. Часто $\alpha = 0.05$, т.е. истинное значение $\theta$ лежит в доверительном интервале с вероятностью $95$%. Часто интервалы строятся не точно, а приближённо. Вместо каких-то параметров подставляются их оценки по выборке. 0. Пусть $X_1, ..., X_n$ - выборка из распределения Бернулли с параметром $p$. $\hat p = \overline{X}$ - оценка параметра $p$. Постройте $(1 - \alpha)$ доверительные интервалы для оценки p:    * используя неравенство Чебышёва.    * асимптотически нормальный доверительный интервал. Решение: $\Exp \hat p = p$ (честно посчитаем мат. ожидание). Неравенство Чебышёва:$$P(\left|\hat p - \Exp \hat p\right| \ge \varepsilon) = P(\left|\hat p - p\right| \ge \varepsilon) \le \frac{\Var \hat p}{\varepsilon^2} = \alpha$$$\alpha$ - вероятность того, что истинное значение не попадёт в д.и. Таким образом:$$\frac{\sqrt{\Var \hat p}}{\sqrt{\alpha}} = \frac{\se \left(\hat p\right)}{\sqrt{\alpha}} = \varepsilon$$Посчитаем дисперсию $\hat p$: $\Var \hat p = \frac{p(1 - p)}{n}$. Поскольку значение $p$ - неизвестно, лучшее, что можно сделать - подставить вместо $p$ его оценку:$$\hat{\Var} \hat p = \frac{\hat p(1 - \hat p)}{n}$$И получить, что $$\hat \varepsilon = \frac{\sqrt{\frac{\hat p(1 - \hat p)}{n}}}{\sqrt{\alpha}}$$Подставив это в неравенство Чебышёва, получим:$$P(\left|\hat p - p\right| \ge \hat \varepsilon) \le \alpha$$$$P(\left|\hat p - p\right| < \hat \varepsilon) \ge 1 - \alpha$$Тогда доверительный интервал выглядит так:$$ \left|\hat p - p\right| < \hat \varepsilon$$$$ - \hat \varepsilon < \hat p - p < \hat \varepsilon$$$$ \hat p - \hat \varepsilon < p < \hat p + \hat \varepsilon$$$$ \hat p - \frac{\sqrt{\frac{\hat p(1 - \hat p)}{n}}}{\sqrt{\alpha}} < p < \hat p + \frac{\sqrt{\frac{\hat p(1 - \hat p)}{n}}}{\sqrt{\alpha}}$$Обратите внимание, д.и. строится именно для $p$, а не для $\hat p$ - мы пытаемся узнать, где находится настоящее значение, а не в какую точку попадёт её оценка по выборке. Теперь построим асимптотически нормальный доверительный интервал. Рассмотрим следующее выражение:$$ \frac{\hat p - p}{\se \left(\hat p\right)}$$Можно показать, что при $n \to +\infty$ эта величина сходится к стандартному нормальному распределению:$$ \frac{\hat p - p}{\se \left(\hat p\right)} \sim N(0, 1)$$Для удобства введём $\xi \sim N(0, 1)$. Тогда $\forall x$:$$P\left(\frac{\hat p - p}{\se \left(\hat p\right)} \le x\right) \approx P\left(\xi \le x\right) = \Phi(x)$$$$ \Phi(x) = \int\limits_{-\infty}^{x}\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}dx$$А теперь построим доверительный интервал:$$ P\left(\left|\frac{\hat p - p}{\se \left(\hat p\right)}\right| \le x\right) \approx P\left(\left|\xi\right| \le x\right) = \Phi(x) - \Phi(-x) = 1 - \alpha$$Т.е. значение лежит в промежутке от $-x$ до $x$ с вероятностью $1 - \alpha$.Этот $x$ находится так:$$ x = \Phi^{-1} \left(1 - \frac{\alpha}{2}\right)$$Теперь можно записать доверительный интервал:$$\left|\frac{\hat p - p}{\se \left(\hat p\right)}\right| \le \Phi^{-1} \left(1 - \frac{\alpha}{2}\right)$$Проделав преобразования, подобно предыдущему пункту задачи, получим:$$\hat p - \se \left(\hat p\right) \cdot \Phi^{-1} \left(1 - \frac{\alpha}{2}\right) \le p \le\hat p + \se \left(\hat p\right) \cdot \Phi^{-1} \left(1 - \frac{\alpha}{2}\right) $$Поскольку $\se \left(\hat p\right)$ зависит от неизвестного $p$, вместо $\se$ подставим $\estSe$:$$\hat p - \sqrt{\frac{\hat p(1 - \hat p)}{n}} \cdot \Phi^{-1} \left(1 - \frac{\alpha}{2}\right) \le p \le\hat p + \sqrt{\frac{\hat p(1 - \hat p)}{n}} \cdot \Phi^{-1} \left(1 - \frac{\alpha}{2}\right) $$Это и есть асимптотически нормальный доверительный интервал. Импортируем нужные модули.<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

%matplotlib inline
import scipy.stats
from statsmodels.distributions.empirical_distribution import ECDF<jupyter_output><empty_output><jupyter_text>1.  Дана выборка из экспоненциального распределения с параметром $\lambda$ (Функция распределения: $F(x) = 1 - e^{-\lambda x}$ при $x > 0$, иначе $0$). $\lambda$ считается неизвестным, для генерации данных используем $\lambda = 3$. Используйте seed = 0.<jupyter_code>lambda_param = 3
random_state = np.random.RandomState(0)<jupyter_output><empty_output><jupyter_text>* Сгенерируйте выборку размера $n=1000$ из экспоненциального распределения.<jupyter_code>n = 1000
sample = random_state.exponential(1 / lambda_param, n)<jupyter_output><empty_output><jupyter_text>* Найдите выборочное среднее и выборочную дисперсию.<jupyter_code>print(f'Выборочное среднее: {np.mean(sample)}')
print(f'Выборочная дисперсия: {np.std(sample) ** 2}')<jupyter_output>Выборочное среднее: 0.33451340292023635
Выборочная дисперсия: 0.11767045932529593<jupyter_text>Постройте графики функции распределения и эмпирической функции распределения при заданном выше $n$. Что можно сказать, глядя на них?<jupyter_code>ecdf = ECDF(sample)
cdf = scipy.stats.expon(scale=1 / lambda_param).cdf(ecdf.x)

plt.figure(figsize = (16, 8))
plt.plot(ecdf.x, ecdf.y, label='Эмпирическая функция распределения', color='blue')
plt.plot(ecdf.x, cdf, color='red', label='Истинная функция распределения')
plt.title('Сравнение истинной и эмпирической функции распределения')
plt.xlabel('x')
plt.ylabel('$\hat F(x)$')
plt.grid()
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>* Постройте $(1-\alpha)$ доверительный интервал для $\lambda$ с помощью оценки $\hat \lambda = \frac{1}{\overline{X}}$.<jupyter_code>alpha = 0.05<jupyter_output><empty_output><jupyter_text>Сперва воспользуемся неравенством Чебышёва: $$P(\left|\tilde \lambda - \Exp \tilde \lambda\right| \ge \varepsilon) \le \frac{\Var \tilde \lambda}{\varepsilon^2} = \alpha$$ Обратите внимание: здесь тильда, а не крышка. В прошлом примере оценка параметра была несмещённой, поэтому в неравенстве Чебышёва мы просто заменили $\Exp \hat p$ на $p$ (они равны). Здесь оценка $\hat \lambda = \frac{1}{\overline{X}}$, приведённая в условии является смещённой. А именно, можно показать, что$$ \Exp \hat \lambda = \Exp \frac{1}{\overline{X}} = \frac{n}{n - 1} \lambda$$Наша цель - найти интервал для $\lambda$, а не для $\hat \lambda$, поэтому можно использовать "исправленную" оценку:$$ \tilde \lambda = \frac{n - 1}{n} \hat \lambda = \frac{n - 1}{n} \frac{1}{\overline{X}}$$ Коэффициент специально выбран так, чтобы оценка была несмещённой. Чтобы решить задачу, нужно найти дисперсию $\tilde \lambda$. Можно показать, что$$\Var \tilde \lambda = \frac{\lambda^2}{n - 2}$$ Таким образом, $\varepsilon = \frac{\lambda}{\sqrt{\alpha \lp n - 2\rp}}$ (смотри предыдущий пример). Теперь чтобы построить доверительный интервал, нужно как-то приблизить $\lambda$. Лучше для этого взять какую-нибудь состоятельную оценку. Взяв $\tilde \lambda$ и воспользовавшись рассуждениями из предыдущего примера, получим:$$ \frac{n - 1}{n} \frac{1}{\overline{X}} - \frac{\frac{n - 1}{n} \frac{1}{\overline{X}}}{\sqrt{\alpha \lp n - 2\rp}} < \lambda < \frac{n - 1}{n} \frac{1}{\overline{X}} + \frac{\frac{n - 1}{n} \frac{1}{\overline{X}}}{\sqrt{\alpha \lp n - 2\rp}}$$ Аналогичным образом получается асимптотически нормальный доверительный интервал: $$ \frac{n - 1}{n} \frac{1}{\overline{X}} - \frac{\frac{n - 1}{n} \frac{1}{\overline{X}}}{\sqrt{n - 2}} \cdot \Phi^{-1} \left(1 -\frac{\alpha}{2}\right) < \lambda < \frac{n - 1}{n} \frac{1}{\overline{X}} + \frac{\frac{n - 1}{n} \frac{1}{\overline{X}}}{\sqrt{n - 2}} \cdot \Phi^{-1} \left(1 -\frac{\alpha}{2}\right) $$ Теперь напишем код, выдающий доверительные интервалы:<jupyter_code>mean = np.mean(sample)
est_lambda = (n - 1) / (n * mean)
est_se = est_lambda / (n - 2) ** 0.5

left_bound = est_lambda - est_se / alpha ** 0.5
right_bound = est_lambda + est_se / alpha ** 0.5
print(f'доверительный интервал на основе неравенства Чебышёва '
      f'для α = {alpha}: ({left_bound}; {right_bound})')

z_alpha = scipy.stats.norm.ppf(1 - alpha / 2, loc=0.0, scale=1.0)
    # 𝛼-квантиль нормального распределения
print(f'𝛼-квантиль нормального распределения: {z_alpha}')
left_bound = est_lambda - z_alpha * est_se
right_bound = est_lambda + z_alpha * est_se
print(f'асимптотически нормальный интервал '
      f'для α = {alpha}: ({left_bound}; {right_bound})')<jupyter_output>доверительный интервал на основе неравенства Чебышёва для α = 0.05: (2.5636598277493; 3.409195019195104)
𝛼-квантиль нормального распределения: 1.959963984540054
асимптотически нормальный интервал для α = 0.05: (2.8011447599339565; 3.1717100870104478)<jupyter_text>Как можно заметить, истинное значение попадает в доверительный интервал. При этом, интервал на основе неравенства Чебышёва несколько шире. * Попадает ли истинное значение в доверительный интервал с ростом $n$? Нужно написать код так, чтобы он работал за $O(n)$ а не $O(n^2)$.<jupyter_code>n_start = 3
ns = np.arange(n_start, n + 1, 1)  # размеры префиксов

means = np.cumsum(sample)[n_start - 1:] / ns
real_values = np.tile(np.array([lambda_param, ]), ns.size)
multiplier = (ns - 1) / ns
stds = multiplier / means / (ns - 2) ** 0.5

estimators = multiplier / means

cheb_lower = estimators - stds / alpha ** 0.5
cheb_upper = estimators + stds / alpha ** 0.5

z_alpha = scipy.stats.norm.ppf(1 - alpha / 2, loc=0.0, scale=1.0)
asymptotic_lower = estimators - stds * z_alpha
asymptotic_upper = estimators + stds * z_alpha

plt.figure(figsize = (16, 8))
plt.plot(ns, real_values, color='red', label='Истинное значение')
plt.plot(ns, cheb_upper, color='green')
plt.plot(ns, cheb_lower, color='green', label='Доверительный интервал с помощью неравенства Чебышёва')

plt.fill_between(ns, asymptotic_upper, asymptotic_lower, 
                 color='blue', label='Асимптотически нормальный доверительный интервад')

plt.plot(ns, estimators, color='orange', label='Несмещённая оценка')

plt.title('Доверительные интервалы')
plt.xlabel('$n$')
plt.ylabel('$\lambda$')
plt.grid()
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Настоящее значение попадает в доверительный интервал. Оценка сходится к истинному значению (что не противоречит её состоятельности).  Самостоятельно/домашка 2. Рассмотрим распределение с плотностью $f(x)$:	$$		f(x) = \begin{cases}			\frac{\beta \alpha^{\beta}}{x^{\beta + 1}}, x \ge \alpha\\			0, x < \alpha		\end{cases}	$$	Здесь $\alpha > 0$ и $\beta > 0$. Пусть $\alpha = 2$ - известно, а $\beta = 5$ - не известно.* (0.5 балла) Сгенерируйте выборку размера $10^4$ из данного распределения. Найдите выборочное среднее и дисперсию.* (0.5 балла) Постройте графики функции распределения и выборочной функции распределения.* (2 балла) Рассмотрим оценку максимального правдоподобия для $\beta$. Используя её, постройте доверительный интервал для $\beta$ на основе неравенства Чебышёва и асимптотически нормальный интервал.* (0.5 балл) Сходится ли оценка максимума правдоподобия к истинному значению $\beta$? 3. Рассмотрим выборку из логнормального распределения с плотностью:$$ f \lp x \rp = \frac{1}{x \sqrt{2 \pi}} e^{-\frac{\lp \ln x - \mu \rp^2}{2}}$$ $\mu = 10$ - неизвестный параметр.* (0.5 балла) Сгенерируйте выборку размера $10^4$ из данного распределения. Найдите выборочное среднее и дисперсию.* (0.5 балла) Постройте графики функции распределения и выборочной функции распределения.* (2 балла) Рассмотрим оценку максимального правдоподобия для $\mu$. Используя её, постройте доверительный интервал для $\mu$ на основе неравенства Чебышёва и асимптотически нормальный интервал.* (0.5 балла) Сходится ли оценка максимума правдоподобия к истинному значению $\mu$? 4. Рассмотрим датасет fifa19. Предположим, что произвольный футболист забивает пенальти с вероятностью $p$. Рассмотрим футболистов из какой-нибудь страны. Теперь сделаем предположение, что $p$ - случайная величина с нормальным распределением с неизвестным математическим ожиданием $\mu$ и неизвестной дисперсией $\sigma^2$. В принципе, $p$ может быть $ 1$, но мы считаем, что нам везёт, и $0 \le p \le 1$ во всех наших наблюдениях. Будем рассматривать футболистов из аргентины.<jupyter_code>data = pd.read_csv('fifa19.zip')
team_name = 'Argentina'
values = data[data['Nationality'] == team_name]['Penalties'].values
only_team = values[~np.isnan(values)]
penalties = only_team / 100<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15789473684210525, 'fraction_numerical': 0.024392120652508464, 'mean_word_length': 4.611830742659758, 'pattern_counts': {'":': 0, '<': 48, '<?xml version=': 0, '>': 38, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '148193', 'license': 'Apache-2.0', 'n_tokens_mistral': 5466, 'n_tokens_neox': 5863, 'n_words': 1332, 'path': 'prac08/stat07_pro.ipynb', 'repo_name': 'tewboat/stat2022', 'size': '12996'}