{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Машинное обучение"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Факультет математики НИУ ВШЭ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2018-2019 учебный год"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Лектор: Илья Щуров\n",
    "\n",
    "Семинаристы: Евгения Ческидова, Евгений Ковалев\n",
    "\n",
    "Ассистенты: Константин Ваниев, Софья Дымченко"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Семинар 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "На этом семинаре мы:\n",
    "\n",
    "- посмотрим, как можно проводить эксперименты из области математической статистики с помощью Python\n",
    "- познакомимся со статистической библиотекой SciPy\n",
    "- проверим работу центральной предельной теоремы\n",
    "- будем проводить собственные исследования!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 Знакомство с SciPy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "С помощью языка Python можно проводить эксперименты из области математической статистики. Например, библиотека для научных и инженерных расчетов [SciPy](http://scipy.github.io/devdocs/index.html) содержит модуль [stats](http://scipy.github.io/devdocs/stats.html#module-scipy.stats), позволяющий работать с распределениями и статистическими функциями. Посмотрим, что можно делать с помощью него."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.stats as sts\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "plt.style.use('ggplot')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# чтобы заработало, почему-то нужно писать в другой ячейке, отдельно от импортирования библиотек\n",
    "plt.rcParams[\"figure.figsize\"] = (10,7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Зададим нормально распределенную случайную величину с параметрами $\\mu = 1$ и $\\sigma = 0.5$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu = 1\n",
    "sigma = 0.5\n",
    "\n",
    "# loc - параметр среднего, scale - параметр среднеквадратичного отклонения\n",
    "norm_rv = sts.norm(loc=mu, scale=sigma)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Построим график функции распределения:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.linspace(-1, 3, 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-1, 3, 100)\n",
    "norm_cdf = norm_rv.cdf(x)\n",
    "#plt.figure(figsize=(10,7))\n",
    "plt.plot(x, norm_cdf)\n",
    "plt.title('Normal RV CDF')\n",
    "plt.xlabel('$x$', fontsize=16)\n",
    "plt.ylabel('$F(x)$', fontsize=16)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Давайте добавим еще распределений!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# равномерное непрерывное\n",
    "a = -1\n",
    "b = 3\n",
    "# loc - параметр левой границы, scale - параметр масштаба (то есть правая граница - loc+scale)\n",
    "uniform_rv = sts.uniform(loc=a, scale=b-a)\n",
    "\n",
    "# биномиальное\n",
    "n = 5\n",
    "p = 0.5\n",
    "binom_rv = sts.binom(n=n, p=p)\n",
    "\n",
    "# Пуассона\n",
    "mu = 5\n",
    "poisson_rv = sts.poisson(mu=mu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-1, 10, 10000)\n",
    "norm_cdf = norm_rv.cdf(x)\n",
    "uniform_cdf = uniform_rv.cdf(x)\n",
    "binom_cdf = binom_rv.cdf(x)\n",
    "poisson_cdf = poisson_rv.cdf(x)\n",
    "plt.plot(x, norm_cdf, label='normal')\n",
    "plt.plot(x, uniform_cdf, label='uniform')\n",
    "plt.plot(x, binom_cdf, label='binomial')\n",
    "plt.plot(x, poisson_cdf, label='poisson')\n",
    "plt.title('Multiple RV CDFs')\n",
    "plt.xlabel('$x$', fontsize=16)\n",
    "plt.ylabel('$F(x)$', fontsize=16)\n",
    "plt.legend(loc='best')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-1, 10, 10000)\n",
    "norm_pdf = norm_rv.pdf(x)\n",
    "uniform_pdf = uniform_rv.pdf(x)\n",
    "binom_pdf = binom_rv.pmf(x)\n",
    "poisson_pdf = poisson_rv.pmf(x)\n",
    "plt.plot(x, norm_pdf, label='normal')\n",
    "plt.plot(x, uniform_pdf, label='uniform')\n",
    "plt.plot(x, binom_pdf, label='binomial')\n",
    "plt.plot(x, poisson_pdf, label='poisson')\n",
    "plt.title('Multiple RV PDFs & PMFs')\n",
    "plt.xlabel('$x$', fontsize=16)\n",
    "plt.ylabel('$p(x)$', fontsize=16)\n",
    "plt.legend(loc='best')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Сгенерируем выборку нормальной случайной величины, упомянутой ранее:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# выбираем число, задающее генератор случайных чисел, чтобы всегда получать один и тот же результат\n",
    "np.random.seed(0)\n",
    "\n",
    "norm_sample = norm_rv.rvs(size=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посчитайте с помощью этой выборки доверительный интервал для среднего с уровнем доверия $95\\%$. Исходя из теоретических соображений, доверительный интервал с уровнем доверия $1 - \\alpha$ может быть представлен в следующем виде:\n",
    "$$\n",
    "\\left(\\hat{\\mu} - z_{1 - \\frac{\\alpha}{2}}\\frac{\\hat{se}\\left(\\hat{\\mu}\\right)}{\\sqrt{n}}, \\ \\hat{\\mu} + z_{1 - \\frac{\\alpha}{2}}\\frac{\\hat{se}\\left(\\hat{\\mu}\\right)}{\\sqrt{n}}\\right)\n",
    "$$\n",
    "Оценку $\\hat{\\mu}$ можно получить методом максимального правдоподобия - она будет равна выборочному среднему. Аналогично, оценка стандартной ошибки $\\hat{se}\\left(\\hat{\\mu}\\right)$ будет выборочным стандартным отклонением. Квантиль же в нашем случае равен $z_{0.975}$.\n",
    "\n",
    "_Квантиль можно посчитать с помощью функции scipy.stats.norm.ppf()_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def CI_calculate(sample):\n",
    "    # YOUR CODE HERE\n",
    "    return (CI_left_bound, CI_right_bound)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Содержит ли построенный интервал искомое значение среднего? А если провести эксперимент много раз, какая доля построенных доверительных интервалов будет содержать это значение? Убедитесь в этом на практике."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# YOUR CODE HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(btw, прикольная визуализация: http://rpsychologist.com/d3/CI/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Проделайте похожий опыт, но с одновыборочным t-критерием Стьюдента. Сгенерируйте выборку той же нормальной случайной величины и проверьте гипотезу о равенстве средних.\n",
    "\n",
    "_Таблица критических значений t-критерия Стьюдента:_ http://www.statisticshowto.com/tables/t-distribution-table/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# YOUR CODE HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "А что будет, если повторить эксперимент много раз? Какова доля (ошибочно) отвергнутых нулевых гипотез?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# YOUR CODE HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2 Проверка ЦПТ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Проверим работу центральной предельной теоремы. Выберем какое-нибудь распределение, не очень похожее на нормальное - например, [Распределение хи-квадрат](https://ru.wikipedia.org/wiki/%D0%A0%D0%B0%D1%81%D0%BF%D1%80%D0%B5%D0%B4%D0%B5%D0%BB%D0%B5%D0%BD%D0%B8%D0%B5_%D1%85%D0%B8-%D0%BA%D0%B2%D0%B0%D0%B4%D1%80%D0%B0%D1%82). Построим случайную величину, имеющую данное распределение, с заданными параметрами:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# число степеней свободы\n",
    "df = 4\n",
    "\n",
    "chi2_rv = sts.chi2(df=df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Сгенерируем выборку размера 1000:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample = chi2_rv.rvs(size=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Построим гистограмму выборки, изобразив поверх нее теоретическую плотность распределения рассматриваемой случайной величины:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# сначала строим гистограмму выборки\n",
    "plt.hist(sample, bins=20, density=True, label='sample points')\n",
    "# теперь строим график плотности\n",
    "x = np.linspace(0, 15, 10000)\n",
    "chi2_pdf = chi2_rv.pdf(x)\n",
    "plt.plot(x, chi2_pdf, label='density')\n",
    "plt.xlabel('$x$')\n",
    "plt.ylabel('number of sample points; $p(x)$')\n",
    "plt.legend(loc='best')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Согласно ЦПТ, выборочное среднее выборки объема $n$ будет иметь распределение, близкое к $N\\left(\\mu, \\frac{\\sigma^2}{n}\\right)$. Мы рассматриваем распределение хи-квадрат с $k=4$ степенями свободы. В этом случае $\\mu = k = 4$, $\\sigma^2 = 2k = 8$. Таким образом, $\\bar{X}_n \\to N\\left(4, \\frac{8}{n}\\right)$ по распределению. \n",
    "\n",
    "Теперь проверим это на практике. Будем генерировать набор выборок размера $n$ (например, 1000 выборок размера $n$), считать набор их выборочных средних, строить их гистограммы и соотносить их с теоретической плотностью соответствующего нормального распределения. Для этого напишем функцию:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample_mean_distr(n):\n",
    "    # YOUR CODE HERE\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим, что получается при подстановке различных значений $n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_mean_distr(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_mean_distr(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_mean_distr(25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_mean_distr(50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_mean_distr(100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Какой вывод можно сделать из этого небольшого проделанного исследования?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 Действительно ли алкоголь влияет на учебу?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "На прошлом семинаре мы сделали предварительное предположение о том, что употребление алкоголя влияет на успехи в учебе: те, кто употребляет больше алкоголя в будние дни, учатся хуже, чем те, кто лучше себя в этом плане контролирует. Проверьте это предположение с помощью t-критерия Стьюдента."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv('math_students.csv', delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# YOUR CODE HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4 Квартет Энскомба"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Рассмотрим четыре различных набора пар $\\left(x_n, y_n\\right)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv('anscombe.csv', index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посчитайте выборочное среднее и выборочную дисперсию для каждого столбца, **не пользуясь функциями .mean(), .std() и .var()** (и проверьте свой ответ с помощью них)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for col in data.columns:\n",
    "    # YOUR CODE HERE\n",
    "    print(col, col_mean, col_var)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Найдите коэффициент корреляции Пирсона для каждой пары $(x_n, y_n)$, __не пользуясь функциями .corr(), .pearsonr() и тому подобными__ (и проверьте свой ответ с помощью них)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(4):\n",
    "    # YOUR CODE HERE\n",
    "    print(x_name, y_name, pearson_corr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Теперь изобразите диаграмму рассеяния для каждой из пар $\\left(x_n, y_n\\right)$ вместе с прямой линейной регрессии $y = 3 + 0.5x$, обозначенной другим цветом."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(4):\n",
    "    # YOUR CODE HERE\n",
    "    plt.title('scatterplot #' + str(i + 1))\n",
    "    plt.xlabel(x_name)\n",
    "    plt.ylabel(y_name)\n",
    "    plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}