{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import pandas as pd\n",
"import scipy.stats as sps\n",
"import matplotlib.pyplot as plt \n",
"import numpy as np\n",
"import seaborn as sns\n",
"\n",
"import random\n",
"# зафиксируем сид для воспроизводимости генерации\n",
"SEED = 42\n",
"np.random.seed(SEED)\n",
"random.seed(SEED)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Критерии вида t-test\n",
"\n",
"### Одновыборочный\n",
"\n",
"Дана одна нормальная выборка $X_1, ..., X_n \\sim \\mathcal{N}(a, \\sigma^2)$.\n",
"\n",
"Критерий проверяет гипотезы\n",
"\n",
"$\\mathsf{H}_0\\colon a = a_0$\n",
"\n",
"$\\mathsf{H}_1\\colon a \\not= a_0$ \n",
"\n",
"`ttest_1samp``(a, popmean): statistic, pvalue`\n",
"\n",
"* `a` — выборка\n",
"* `popmean` — равно $a_0$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Сгенерируйте выборку $X_1, ..., X_n \\sim \\mathcal{N}(0, 1)$ и с помощью критерия проверьте:\n",
"- равенство среднего нулю \n",
"- равенство среднего 0.5"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"size=100\n",
"\n",
"sample = <...>\n",
"<...>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Двухвыборочный\n",
"\n",
"#### Независимые выборки\n",
"\n",
"Даны две независимые нормальные выборки\n",
"\n",
"* $X_1, ..., X_n \\sim \\mathcal{N}(a_1, \\sigma_1^2)$,\n",
"\n",
"* $Y_1, ..., Y_m \\sim \\mathcal{N}(a_2, \\sigma_2^2)$.\n",
"\n",
"Критерий проверяет для их гипотезы о равенстве среднего:\n",
"\n",
"$\\mathsf{H}_0\\colon a_1 = a_2$\n",
"\n",
"$\\mathsf{H}_1\\colon a_1 \\not= a_2$ \n",
"\n",
"`ttest_ind``(a, b, equal_var=True): statistic, pvalue`\n",
"\n",
"`a`, `b` — выборка\n",
"\n",
"`equal_var` — известно ли равенство дисперсий"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Сгенерируйте выборки $X_1, ..., X_n \\sim \\mathcal{N}(0, 1)$ и $X_1, ..., X_m \\sim \\mathcal{N}(1, 1)$. Используя критерий, проверьте равенство средних двух выборок."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sample_1 = <...>\n",
"sample_2 = <...>\n",
"<...>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Сгенерируйте выборки $X_1, ..., X_n \\sim \\mathcal{N}(0, 1)$ и $X_1, ..., X_m \\sim \\mathcal{N}(1, 7)$. Используя критерий, проверьте равенство средних двух выборок."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sample_1 = <...>\n",
"sample_2 = <...>\n",
"<...>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Связные выборки\n",
"\n",
"Даны две связные нормальные выборки\n",
"\n",
"* $X_1, ..., X_n \\sim \\mathcal{N}(a_1, \\sigma_1^2)$,\n",
"\n",
"* $Y_1, ..., Y_n \\sim \\mathcal{N}(a_2, \\sigma_2^2)$.\n",
"\n",
"\n",
"Критерий проверяет для их гипотезы о равенстве среднего:\n",
"\n",
"$\\mathsf{H}_0\\colon a_1 = a_2$\n",
"\n",
"$\\mathsf{H}_1\\colon a_1 \\not= a_2$ \n",
"\n",
"`ttest_rel``(a, b): statistic, pvalue`\n",
"\n",
"`a`, `b` — выборка"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Сгенерируйте выборку $X_1, ..., X_n \\sim \\mathcal{N}(0, 1)$. Вторую сгенерируйте по формуле выборка_1 + случайный шум, случайный шум из $\\mathcal{N}(0, 0.5)$. Используя критерий, проверьте гипотезу о равенстве среднего."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sample_1 = <...>\n",
"sample_2 = <...>\n",
"<...>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Сгенерируйте выборку $X_1, ..., X_n \\sim \\mathcal{N}(0, 1)$. Вторую сгенерируйте по формуле выборка_1 + случайный шум, случайный шум из $\\mathcal{N}(0.5, 0.5)$. Используя критерий, проверьте гипотезу о равенстве среднего."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sample_1 = <...>\n",
"sample_2 = <...>\n",
"<...>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Пример: ирисы Фишера\n",
"\n",
"Визуализация данных"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"df = sns.load_dataset(\"iris\")\n",
"\n",
"g = sns.PairGrid(df, hue='species')\n",
"g.map_lower(sns.kdeplot, cmap =\"Blues_d\")\n",
"g.map_upper(plt.scatter)\n",
"g.map_diag(sns.kdeplot, lw=3);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Как выглядят данные"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"df.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Виды ирисов"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"np.unique(df.species)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sps.ttest_ind(df[df.species == 'setosa'].sepal_length, \n",
" df[df.species == 'versicolor'].sepal_length,\n",
" equal_var=False)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sps.ttest_ind(df[df.species == 'virginica'].sepal_length, \n",
" df[df.species == 'versicolor'].sepal_length,\n",
" equal_var=False)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sps.ttest_ind(df[df.species == 'virginica'].sepal_width, \n",
" df[df.species == 'versicolor'].sepal_width,\n",
" equal_var=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*Замечание.* Строго говоря, неоходима поправка на множественное тестирование гипотез.\n",
"\n",
"**Вывод**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Множественная проверка гипотез"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"С помощью статистических методов можно проверить человека на наличие экстрасенсорных способностей: предложим ему угадать последовательность, состоящую из двух цветов, длины 100. \n",
"\n",
"Сформулируем задачу на статистическом языке: \n",
"\n",
"$X_1...X_{100}$ — выборка из распределения $Bern(p)$\n",
"\n",
"$p=0.5$ отвечает случайному угадыванию.\n",
"\n",
"Проверьте гипотезу: $\\mathsf{H}_0 \\colon p=0.5$ vs $\\mathsf{H}_1 \\colon p \\neq 0.5$. Используйте критерий Вальда.\n",
"\n",
"В качестве асимптотически нормальной оценки можно использовать $\\widehat{p} = \\overline{X}$ с асимптотической дисперсией $\\sigma^2(p) = p (1 - p)$.\n",
"\n",
"Выпишем состоятельную оценку дисперсии и статистику критерия Вальда:\n",
"\n",
"$\\widehat{\\sigma} = \\sqrt{\\overline{X} (1 - \\overline{X})}$, $W = \\sqrt{n} \\frac{\\overline{X} - 0.5}{\\sqrt{\\overline{X} (1 - \\overline{X})}}$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Оценим реальный уровень значимости для этого критерия при размере выборки равном 100. К чему он должен быть близок? Для скорости вычислений используйте количество выборок равное $10^3$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sample_size = 100\n",
"sample_count = 1000\n",
"\n",
"theta = 0.5"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def wald_test(sample, theta, estimation_theta, estimation_sigma, alternative='two_sided'):\n",
" \"\"\"\n",
" param sample: реализация выборки\n",
" param theta: истинное значение параметра\n",
" param estimation_theta: оценка параметра\n",
" param estimation_sigma: оценка асимптотической дисперсии оценки estimation_sigma\n",
" param alternative: вид альтернативной гипотезы, может принимать одно из значений 'two_sided', 'less', 'greater'\n",
"\n",
" return statistic\n",
" return p_value\n",
" \"\"\"\n",
"\n",
" alpha = 0.05\n",
" z = sps.norm.ppf(1 - alpha/2)\n",
" n = len(sample)\n",
" statistic = np.sqrt(n) * (estimation_theta - theta) / estimation_sigma\n",
"\n",
" if alternative == 'two_sided':\n",
" p_value = sps.norm.sf(np.abs(statistic)) + sps.norm.cdf(-np.abs(statistic))\n",
" conf_int = round(estimation_theta - z*estimation_sigma/np.sqrt(n), 4), round(estimation_theta + z*estimation_sigma/np.sqrt(n), 4)\n",
"\n",
"\n",
" elif alternative == 'less':\n",
" p_value = sps.norm.cdf(statistic)\n",
" conf_int = (-np.inf, round(estimation_theta + z*estimation_sigma/np.sqrt(n), 4))\n",
"\n",
" \n",
" elif alternative == 'greater':\n",
" p_value = sps.norm.sf(statistic)\n",
" conf_int = (round(estimation_theta + z*estimation_sigma/np.sqrt(n), 4), np.inf)\n",
"\n",
" else:\n",
" raise ValueError('alternative name is wrong')\n",
"\n",
" return statistic, p_value, conf_int"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Оценим реальный уровень значимости"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sample = <...>\n",
"\n",
"estimation_theta = <...>\n",
"estimation_sigma = <...>\n",
"\n",
"\n",
"counter = 0\n",
"\n",
"for i in range(sample_count):\n",
" _, p_value, conf_int = wald_test(sample[i], theta, estimation_theta[i], estimation_sigma[i])\n",
" is_rejected = <...>\n",
" if is_rejected:\n",
" counter += 1\n",
" \n",
"counter / sample_count\n",
" \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Теперь представим, что мы хотим проверить большое количество людей на экстрасенсорные способности с помощью данного критерия.\n",
"\n",
"Проведите аналогичный эксперимент: сгенерируйте $10^3$ выборок размера $100$ для $100$ людей. Посчитайте, сколько раз из 1000 в вашем наборе из 100 выборок хотя бы для одной гипотеза будет отвергнута. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sample_all = <...>\n",
"\n",
"counter = 0\n",
"for sample in sample_all:\n",
" estimation_theta = <...>\n",
" estimation_sigma = <...>\n",
"\n",
" for i in range(100):\n",
" _, p_value, conf_int = wald_test(sample[i], theta, estimation_theta[i], estimation_sigma[i])\n",
" is_rejected = <...>\n",
" if is_rejected:\n",
" counter += 1\n",
" break"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"counter / 1000"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Вывод:** <...>\n",
"\n",
"На лекции вы прошли метод, позволяющий не накапливать ошибку 1 рода. В этом методе необходимо использовать уровень значимости, зависящий от количества проверяемых одновременно гипотез. \n",
"\n",
"**Чему равен этот уровень значимости, если одновременно проверяются n гипотез?**\n",
"\n",
"**Ответ:** <...>\n",
"\n",
"Проведите предыдущий эксперимент с использованием корректной процедуры. Поскольку в реализованной выше функции $\\alpha$ зафиксировано, используйте критерий отвержения гипотезы с помощью p-value."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sample_all = <...>\n",
"\n",
"counter = 0\n",
"for sample in sample_all:\n",
" estimation_theta = <...>\n",
" estimation_sigma = <...>\n",
"\n",
" for i in range(100):\n",
" _, p_value, conf_int = wald_test(sample[i], theta, estimation_theta[i], estimation_sigma[i])\n",
" is_rejected = <...>\n",
" if is_rejected:\n",
" counter += 1\n",
" break "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"counter / 1000"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Вывод:** <...>"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.12"
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 2
}