{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import scipy.stats as sst\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "В scipy.stats используются следующие функции:\n", "\n", "Пусть $X$ — класс, реализующий некоторое распределение. Конкретное распределение с параметрами `params` можно получить как **`X(params)`**. \n", "У него доступны следующие методы:\n", "\n", "* `X(params).rvs(size=N)` — генерация выборки размера $N$ (Random VariateS). Возвращает `numpy.array`;\n", "* `X(params).cdf(x)` — значение функции распределения в точке $x$ (Cumulative Distribution Function);\n", "* `X(params).ppf(q)` — $q$-квантиль (Percent Point Function);\n", "* `X(params).mean()` — математическое ожидание;\n", "* `X(params).median()` — медиана ($1/2$-квантиль);\n", "* `X(params).var()` — дисперсия (Variance);\n", "* `X(params).std()` — стандартное отклонение = корень из дисперсии (Standard Deviation).\n", "* `X(params).pdf(x)` — для непрерывных распределений значение плотности в точке $x$ (Probability Density Function);\n", "* `X(params).pmf(k)` — для дискретных распределений значение дискретной плотности в точке $k$ (Probability Mass Function);\n", "\n", "Эти методы применимы как к конкретному распределению `X(params)`, так и к самому классу `X`. Во втором случае параметры передаются в сам метод. Например, вызов `X.rvs(size=N, params)` эквивалентен `X(params).rvs(size=N)`. При работе с распределениями и случайными величинами рекомендуем использовать первый способ, посколько он больше согласуется с математическим синтаксисом теории вероятностей.\n", "\n", "Параметры могут быть следующими:\n", "* `loc` — параметр сдвига;\n", "* `scale` — параметр масштаба;\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# I. Параметрические тесты\n", "\n", "Воспользуемся параметрическим критерием Стьюдента для проверки равенства средних.\n", "\n", "Сгенерируем 3 выборки с разными параметрами нормального распределения - среднее=[1; 3; 1], дисперсия=[1; 1; 5] и размером=[50; 40; 60].\n", "\n", "Для этого воспользуемся вероятностными функциями из numpy" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#\n", "mean = [1, 2, 1]\n", "variance = [1, 1, 4]\n", "size = [50, 40, 60]\n", "\n", "if len(mean) != len(variance) or len(variance) != len(size):\n", " print(\"Ошибка в исходных данных! Размер массивов должен совпадать!\")\n", "#\n", "task_arr = []\n", "for i in range(len(size)):\n", " task_arr.append(np.random.normal(loc=mean[i], scale=np.sqrt(variance[i]), size=size[i]))\n", "#\n", "full_arr = np.array([])\n", "for arr in task_arr:\n", " for num in arr:\n", " full_arr = np.append(full_arr, num)\n", "\n", "print(\"mean of full =\", np.mean(full_arr))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " 1) Рассчитаем t-критерий для средних значений двух независимых выборок - ttest_ind.\n", "\n", "Если equal_var=True, то выполняется стандартный независимый тест с двумя выборками, который предполагает равные дисперсии совокупности. \n", "Если equal_var=False, то выполнится выриант t-критерия Уэлча, который не предполагает равную дисперсию." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "alpha = 0.05\n", "# Вычисляем t-статистику и p-value\n", "t_stat, p_value = sst.ttest_ind(task_arr[0], task_arr[1], equal_var=True)\n", "print(\"t-test -equal variances\")\n", "if p_value < alpha:\n", " print(\"Отвергаем нулевую гипотезу, т.к. существует статистически значимое различие между выборками.\")\n", "else:\n", " print(\"Не отвергаем нулевую гипотезу - не существует статистически значимого различия между выборками.\")\n", "\n", "# Вычисляем t-статистику и p-value\n", "t_stat, p_value = sst.ttest_ind(task_arr[0], task_arr[1], equal_var=False)\n", "print(\"t-test -not equal variances\")\n", "if p_value < alpha:\n", " print(\"Отвергаем нулевую гипотезу - существует статистически значимое различие между выборками.\")\n", "else:\n", " print(\"Не отвергаем нулевую гипотезу - не существует статистически значимого различия между выборками.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Получите результаты сравнения двух других пар выборок - task_arr[0] и task_arr[2]; task_arr[1] и task_arr[2]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Вычисляем t-статистику и p-value" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "2) Используем F-тест для проверки равенства дисперсий.\n", "Вычислим F-статистику и p-value для пары выборок." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Вычисляем степени свободы выборок\n", "x = task_arr[0]\n", "df1 = len(x) - 1\n", "y = task_arr[1]\n", "df2 = len(y) - 1\n", "\n", "# Вычисляем F-статистику и p-value\n", "p_value = sst.f(df1, df2).cdf( x.var() / y.var())\n", "\n", "# Set the significance level\n", "alpha = 0.05\n", "# Сравниваем F-value с уровнем значимости\n", "if p_value < alpha:\n", " print(\"Отвергаем нулевую гипотезу - существует статистически значимое различие в дисперсиях.\")\n", "else:\n", " print(\"Не отвергаем нулевую гипотезу - не существует статистически значимого различия в дисперсиях.\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Получите результаты сравнения двух других пар выборок - task_arr[0] и task_arr[2]; task_arr[1] и task_arr[2]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Вычисляем F-статистику и p-value" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "3) В фреймворке scipy.stats для проверки нормальности распределения также реализован критерий д'Агостиньо-Пирсона." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Генеральная совокупность\n", "#population = np.random.normal(0, 1, size=1000)\n", "\n", "# Применяем критерий Агостиньо-Пирсона для проверки нормальности распределения\n", "statistic, p_value = sst.normaltest(full_arr)\n", "print(\"критерий Агостиньо-Пирсона\")\n", "# Сравниваем p-value с уровнем значимости alpha\n", "if p_value < alpha:\n", " print(\"Отвергаем нулевую гипотезу, генеральная совокупность не распределена нормально.\")\n", "else:\n", " print(\"Не отвергаем нулевую гипотезу, генеральная совокупность распределена нормально.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### II. Метод получения выборки аппроксимацией\n", "\n", "Допустим, у нас нет хорошего простого аналитического распределения, и мы не можем использовать его для аппроксимации нашего распределения. Можно использовать аппроксимацию нашей выборки на базе равномерного распределения методом отбрасывания некоторых точек. Выполняется следующим образом:\n", "\n", "1. Выберем равномерно распределенное 'x'.\n", "2. Выберем равномерно распределенное `y` от `0` до максимальной вероятности в нашем распределении плотности.\n", "3. Если $y > p(x)$ то вычеркнем такую точку.\n", "\n", "Попробуем выполнить выборку из *ненормализованного* распределения $p(x) = \\sin(x^2) + 1$ от $0 \\rightarrow 4$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# сгенерировать выборку из распределения, которого нет в scipy.stats\n", "def my_pdf(x): return np.sin(x**2) + 1\n", "\n", "# Рисуем график плотности непрерывного распределения\n", "# функция linspace - Создание массива из равномерно распределенных значений, \n", "# это последовательности чисел с постоянным шагом, или равномерная сетка в диапазоне значений\n", "xs = np.linspace(0, 4, 200)\n", "ps = my_pdf(xs)\n", "plt.plot(xs, ps)\n", "# закрасим область под графиком\n", "plt.fill_between(xs, 0, ps, alpha=0.1)\n", "plt.xlabel('Значение', fontsize=12)\n", "plt.ylabel('Плотность', fontsize=12)\n", "plt.grid(ls=':')\n", "# ограничим оси\n", "plt.xlim(0, 4); plt.ylim(0, 2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n = 100\n", "# создадим поле равномерно распределенных значений\n", "random_x = sst.uniform.rvs(loc=0, scale=4, size=n)\n", "random_y = sst.uniform.rvs(loc=0, scale=2, size=n)\n", "\n", "plt.scatter(random_x, random_y)\n", "plt.plot(xs, ps, c=\"k\")\n", "# закрасим область под графиком\n", "plt.fill_between(xs, 0, ps, color=\"y\", alpha=0.1)\n", "plt.xlabel('Значение', fontsize=12)\n", "plt.ylabel('Плотность', fontsize=12)\n", "plt.grid(ls=':')\n", "# ограничим оси\n", "plt.xlim(0, 4), plt.ylim(0, 2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# отфильтруем точки по нашей функции\n", "passed = random_y <= my_pdf(random_x)\n", "# отобразим две группы точек\n", "plt.scatter(random_x[passed], random_y[passed])\n", "plt.scatter(random_x[~passed], random_y[~passed], marker=\"x\", s=30, alpha=0.5)\n", "\n", "plt.plot(xs, ps, c=\"w\")\n", "# закрасим область под графиком\n", "plt.fill_between(xs, 0, ps, color=\"k\", alpha=0.1)\n", "plt.xlabel('Значение', fontsize=12)\n", "plt.ylabel('Плотность', fontsize=12)\n", "plt.grid(ls=':')\n", "# ограничим оси\n", "plt.xlim(0, 4), plt.ylim(0, 2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n2 = 100000\n", "x_test = sst.uniform.rvs(scale=4, size=n2)\n", "x_final = x_test[sst.uniform.rvs(scale=2, size=n2) <= my_pdf(x_test)]\n", "print(len(x_final))\n", "\n", "# используем метод интегрирования с помощью выборок вдоль заданной оси и правила Симпсона.\n", "from scipy.integrate import simpson\n", "\n", "plt.hist(x_final, density=True, histtype=\"step\", label=\"гистограмма выборки\")\n", "plt.plot(xs, ps / simpson(ps, x=xs), c=\"k\", label=\"функция плотности выборки\")\n", "plt.xlabel('Значение', fontsize=12)\n", "plt.ylabel('Плотность', fontsize=12)\n", "plt.grid(ls=':')\n", "plt.legend(loc=3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### III. Выборка методом обратного преобразования\n", "\n", "Концептуально так: мы знаем, что для всех PDF значение CDF изменяется от 0 до 1. \n", "Если можно равномерно отобрать CDF от 0 до 1, можем ли мы инвертировать нашу функцию, чтобы восстановить значение $x$, которое дает выборочное значение CDF ? \n", "\n", "Для некоторых функций - это возможно, но для некоторых - нет, математика неразрешима!\n", "\n", "Приведенная выше функция сложна, поэтому найдем функцию, которую проще использовать для сэмплирования.\n", "\n", "Допустим, наш PDF равен $p(x) = 3 x ^ 2$ из $0 \\rightarrow 1$. В отличие от предыдущей функции, эта нормализована.\n", "\n", "Можно найти CDF с помощью\n", "\n", "$$ CDF(x) = \\int_0^x p(x ^\\prime) dx ^\\prime = x ^ 3 $$\n", "\n", "Теперь есть CDF и мы хотим его инвертировать, то есть на данный момент у нас есть значение $x$ и получаем значение $y$ - CDF. \n", "Мы хотим присвоить ему значение CDF и получить значение $x$:\n", "\n", "$$ y = x^3 \\rightarrow x = y^3 $$\n", "\n", "здесь получится $y = x^{1/3}$ или в виде CDF $x = CDF^{1/3}$\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def pdf(x): return 3 * x**2\n", "def cdf(x): return x**3\n", "def icdf(cdf): return cdf**(1/3)\n", "\n", "# Создание массива из равномерно распределенных значений\n", "xs = np.linspace(0, 1, 100)\n", "pdfs = pdf(xs)\n", "cdfs = cdf(xs)\n", "n = 2000\n", "u_samps = sst.uniform.rvs(size=n)\n", "x_samps = icdf(u_samps)\n", "\n", "# построим две диаграммы рядом PDF и CDF\n", "fig, axes = plt.subplots(ncols=2, figsize=(10, 4))\n", "axes[0].plot(xs, pdfs, color=\"k\", label=\"PDF-функция плотности\")\n", "axes[0].hist(x_samps, density=True, histtype=\"step\", label=\"моя выборка\", bins=50)\n", "axes[1].plot(xs, cdfs, color=\"k\", label=\"CDF-функция распределения\")\n", "\n", "axes[1].hlines(u_samps, 0, x_samps, linewidth=0.1, alpha=0.3)\n", "axes[1].vlines(x_samps, 0, u_samps, linewidth=0.1, alpha=0.3)\n", "\n", "axes[0].legend(), axes[1].legend()\n", "axes[0].set_xlim(0, 1), axes[0].set_ylim(0, 3)\n", "axes[1].set_xlim(0, 1), axes[1].set_ylim(0, 1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Проверим возможность применения метода обратного преобразования для получения своего распределения" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# используем интерполяцию из фреймворка scipy\n", "from scipy.interpolate import interp1d\n", "\n", "def my_pdf(x): return np.sin(x**2) + 1\n", "# Создание массива из равномерно распределенных значений\n", "xs = np.linspace(0, 4, 4000)\n", "pdfs = my_pdf(xs)\n", "cdfs = pdfs.cumsum() / pdfs.sum() # опасное действие! результат не гарантирован!\n", "# при получении ошибки повторно выполняйте код до получения результата!\n", "\n", "# Создание выборки из равномерно распределенных значений\n", "u_samps = sst.uniform.rvs(size=4000)\n", "x_samps = interp1d(cdfs, xs)(u_samps)\n", "\n", "# построим две диаграммы рядом PDF и CDF\n", "fig, axes = plt.subplots(ncols=2, figsize=(10,4))\n", "axes[0].hist(x_samps, density=True, histtype=\"step\", label=\"моя гистограмма\", bins=50)\n", "axes[0].plot(xs, pdfs/4.747, color=\"k\", label=\"теор. PDF\")\n", "axes[0].legend(loc=3), axes[0].set_xlim(0, 4)\n", "\n", "axes[1].plot(xs, cdfs, color=\"k\", label=\"эмпирич.CDF\")\n", "axes[1].hlines(u_samps, 0, x_samps, linewidth=0.1, alpha=0.3)\n", "axes[1].vlines(x_samps, 0, u_samps, linewidth=0.1, alpha=0.3)\n", "axes[1].legend(loc=2), axes[1].set_xlim(0, 4), axes[1].set_ylim(0, 1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## IV. Интерполяция данных\n", "\n", "Нам нужен способ интерполяции между точками с некоторой точностью.\n", "Кратко рассмотрим интерполяцию." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Вариант вычисления эмпирических PDF, CDF и SF.\n", "\n", "Предположим, что у нас есть несколько значений из некоторого распределения. Как можно получить непрерывную функцию для дальнейшего использования?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xs = [0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, \n", " 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0]\n", "ys = [0.2, 0.165, 0.167, 0.166, 0.154, 0.134, 0.117, \n", " 0.108, 0.092, 0.06, 0.031, 0.028, 0.048, 0.077, \n", " 0.103, 0.119, 0.119, 0.103, 0.074, 0.038, 0.003]\n", "\n", "plt.scatter(xs, ys)\n", "plt.xlabel(\"x\")\n", "plt.ylabel(\"Опытная (Observed) PDF\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# используем интерполяцию из фреймворка scipy\n", "from scipy.interpolate import interp1d\n", "\n", "x = np.linspace(min(xs), max(xs), 1000)\n", "y1 = interp1d(xs, ys)(x)\n", "y2 = interp1d(xs, ys, kind=\"nearest\")(x)\n", "y3 = interp1d(xs, ys, kind=\"quadratic\")(x)\n", "y4 = interp1d(xs, ys, kind=\"cubic\")(x)\n", "\n", "plt.scatter(xs, ys, s=30, label=\"данные\", c=\"w\")\n", "plt.plot(x, y1, label=\"линейная (по умолчанию)\")\n", "plt.plot(x, y2, label=\"Nearest\", alpha=0.2)\n", "plt.plot(x, y3, label=\"Quadratic\", ls='-')\n", "plt.plot(x, y4, label=\"Cubic\", ls='-')\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Используя `interp1d` можно найти вероятность для любого значения `x`.\n", "\n", "Как можно вычислить CDF и вероятность для значений между границами? - С помощью `scipy.integrate`!\n", "Есть несколько методов вычисления интеграла:\n", "\n", "* `scipy.integrate.trapz` для низкой точности, но высокой скорости вычисления\n", "* `scipy.integrate.simpson` для средней точности, но достаточно высокой скорости вычисления\n", "* `scipy.integrate.quad` для высокой точности, но низкой скорости вычисления\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy.integrate import simpson\n", "\n", "def get_prob(xs, ys, a, b, resolution=1000):\n", " x_norm = np.linspace(min(xs), max(xs), resolution)\n", " y_norm = interp1d(xs, ys, kind=\"quadratic\")(x_norm)\n", " normalisation = simpson(y_norm, x=x_norm)\n", " x_vals = np.linspace(a, b, resolution)\n", " y_vals = interp1d(xs, ys, kind=\"quadratic\")(x_vals)\n", " return simpson(y_vals, x=x_vals) / normalisation\n", "\n", "def get_cdf(xs, ys, v):\n", " return get_prob(xs, ys, min(xs), v)\n", "\n", "# в интервале от 0 до 10\n", "print(get_prob(xs, ys, 0, 10))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# в интервале от 6 до 9.3\n", "v1, v2 = 6, 9.3\n", "area = get_prob(xs, ys, v1, v2)\n", "\n", "plt.scatter(xs, ys, s=30, label=\"данные\", color=\"w\")\n", "plt.plot(x, y3, linestyle=\"-\", label=\"Interpolation\")\n", "plt.fill_between(x, 0, y3, where=(x>=v1)&(x<=v2), alpha=0.2)\n", "plt.annotate(f\"p = {area:.3f}\", (7, 0.05))\n", "plt.legend()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x_new = np.linspace(min(xs), max(xs), 100)\n", "cdf_new = [get_cdf(xs, ys, i) for i in x_new]\n", "cp_cdf = y3.cumsum() / y3.sum()\n", "\n", "plt.plot(x_new, cdf_new, label=\"интерполяция CDF\")\n", "plt.plot(x, cp_cdf, label=\"простая CDF для типовых случаев\")\n", "plt.ylabel(\"CDF\")\n", "plt.xlabel(\"x\")\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Обратите внимание, что метод CDF, который здесь реализован, является простым, но довольно медленным.\n", "\n", "Обычно следует один раз определить функцию CDF с высокой точностью и сохранить результат, а затем можно просто интерполировать его. Пересчитывать его каждый раз - пустая трата процессорного времени. \n", "\n", "Для подавляющего большинства приложений метод `cumsum()` подходит, если у вас есть массив с достаточным количеством точек." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# пример - произвольная функция некоторого распределния\n", "def fn(xs):\n", " return np.exp((xs + 10.5)**0.1) + np.sin(xs) / (xs + 1) + np.cos(2.5 * xs**0.5)**2\n", "\n", "xs = np.array([1, 2, 4, 5, 7, 9, 10])\n", "ys = fn(xs)\n", "# Создание массива из равномерно распределенных значений\n", "x_fine = np.linspace(1, 10, 200)\n", "y_fine = fn(x_fine)\n", "\n", "plt.plot(xs, ys, 'o', label=\"значения\", markersize=5)\n", "plt.plot(x_fine, y_fine, label=\"функция\", color=\"k\", alpha=0.5, linestyle=\"-\")\n", "plt.legend()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# линейная интерполяция\n", "y_linear = interp1d(xs, ys, kind=\"linear\")(x_fine)\n", "# квадратичная интерполяция\n", "y_quad = interp1d(xs, ys, kind=\"quadratic\")(x_fine)\n", "# кубическая интерполяция\n", "y_cubic = interp1d(xs, ys, kind=\"cubic\")(x_fine)\n", "\n", "plt.plot(xs, ys, 'o', label=\"значения\", markersize=5)\n", "plt.plot(x_fine, y_fine, label=\"функция\", color=\"b\", ls=\"-\", alpha=0.75)\n", "plt.plot(x_fine, y_linear, label=\"линейная\")\n", "plt.plot(x_fine, y_quad, label=\"квадрат.\")\n", "plt.plot(x_fine, y_cubic, label=\"кубич.\")\n", "plt.legend()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "class empiric_gen(sst.rv_continuous):\n", " def _pdf(self, x):\n", " return np.exp((x + 10.5)**0.1) + np.sin(x) / (x + 1) + np.cos(2.5 * x**0.5)**2\n", "\n", "empir = empiric_gen(a=1, b=9, name='experiment')\n", "sampl_emp = empir.rvs(size=200)\n", "\n", "print('Первые 10 значений выборки:\\n', sampl_emp[:10])\n", "print(f'Выборочное среднее: {sampl_emp.mean():.3f}' )\n", "print(f'Выборочная дисперсия: {sampl_emp.var():.3f}' )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Если дискретная случайная величина может принимать небольшое число значений, то можно и не создавать новый класс, а явно указать эти значения и их вероятности (как в табличном формате Excel). " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "some_distribution = sst.rv_discrete(\n", " name='мое распределение',\n", " values=([1, 3, 7], [0.6, 0.1, 0.3]) # значения и вероятности\n", ")\n", "sample = some_distribution.rvs(size=200)\n", "\n", "print('Первые 20 значений выборки:\\n', sample[:20])\n", "print(f'Выборочное среднее: {sample.mean():.3f}' )\n", "print('Частота значений по выборке:', (sample == 1).mean(), (sample == 3).mean(), (sample == 7).mean())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# V. Интервальная оценка выборочных данных\n", "\n", "Смоделируем НЕнормальное распределение, чтобы проиллюстрировать центральную предельную теорему (ЦПТ).\n", "\n", "ЦПТ утверждает, что распределение, которое получается из суммы многих выборок произвольного распределения, приближается к нормальному распределению, поскольку суммируется все больше и больше различных точек.\n", "\n", " Функция `my_get_data` реализует способ получения выборки произвольного распределения." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def my_get_data(n):\n", " # объединим данные из двух различных распределений\n", " data = np.concatenate(\n", " (sst.expon.rvs(scale=1, size=n//2),\n", " sst.skewnorm.rvs(5, loc=3, size=n//2)))\n", " # перемешаем значения в массиве\n", " np.random.shuffle(data)\n", " return data\n", "n=2000\n", "plt.hist(my_get_data(n), bins=50)\n", "print(n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Это определенно не нормальное распределение.\n", "\n", "Обычно в практических экспериментах нет возможности запрашивать бесконечные объемы данных - приходится постепенно их собирать, что требует времени, усилий или денег. \n", "\n", "Предположим, получено только 10 точек данных. Вычислим среднее значение этой выборки." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "d10 = my_get_data(10)\n", "print(d10.mean())\n", "# построим гистограмму по 10 выборкам с 10 значениями\n", "ave_means = [my_get_data(10).mean() for i in range(10)]\n", "plt.hist(ave_means, bins=5)\n", "# вычислим доверительный интервал\n", "ti10=sst.t.interval(0.95, loc=np.mean(ave_means), df=len(ave_means)-1, scale=sst.sem(ave_means))\n", "## нарисуем верт.линии по значению среднего и довер.интервал\n", "plt.axvline(np.mean(ave_means), label=\"mean\",c='r', ls=\"--\")\n", "plt.axvline(ti10[0], label=\"mean-\",c='y', ls=\"--\", alpha=0.5)\n", "plt.axvline(ti10[1], label=\"mean+\",c='y', ls=\"--\", alpha=0.5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Если повторить эту операцию несколько раз, то увидим, что ответ меняется после каждого выполнения кода!\n", "\n", "Дело в том, что величина среднего выборочного не случайна. Между прогонами это так, но у этого результата есть определенные статистические свойства. \n", "\n", "Давайте посмотрим, что произойдет, если нам удастся получить 10 точек данных 1000 раз. Проанализируем 1000 результатов mean()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ave_means = [my_get_data(10).mean() for i in range(1000)]\n", "plt.hist(ave_means, bins=50)\n", "print(np.std(ave_means))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Распределение становится все больше и больше похожим на нормальное распределение. Пик распределения приближается к истинному среднему по совокупности, и ширина этого распределения является функцией как ширины совокупности, так и количества использованных нами выборок. \n", "\n", "Если использовать 100 точек данных, можно ожидать меньшего разброса между выборками. Давайте протестируем это предположение." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ave_means = [my_get_data(100).mean() for i in range(1000)]\n", "plt.hist(ave_means, bins=50)\n", "print(np.std(ave_means))\n", "ti10=sst.t.interval(0.95, loc=np.mean(ave_means), df=len(ave_means)-1, scale=sst.sem(ave_means)) \n", "## нарисуем верт.линию по значению среднего\n", "plt.axvline(np.mean(ave_means), label=\"mean\",c='r', ls=\"--\")\n", "plt.axvline(ti10[0], label=\"mean-\",c='y', ls=\"--\", alpha=0.5)\n", "plt.axvline(ti10[1], label=\"mean+\",c='y', ls=\"--\", alpha=0.5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Это распределение выглядит более гауссовым! Можно заметить что и дисперсия здесь меньше, чем ранее." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "num_samps = [10, 50, 100, 500, 1000, 5000, 10000]\n", "my_stda = []\n", "# цикл получения новых данных\n", "for n in num_samps:\n", " my_stda.append(np.std([my_get_data(n).mean() for i in range(1000)]))\n", "# построение точечного графика с меткой \"о\"\n", "plt.plot(num_samps, my_stda, 'o', label=\"значения станд.откл.\")\n", "# построение линейного графика аппроксимации зависимости std\n", "plt.plot(num_samps, 1 / np.sqrt(num_samps), label=\"аппрокс.функция\", alpha=0.5)\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Случайная функция не произвольна, она равна $1/\\sqrt{N}$.\n", "\n", "Распределение средних имеет стандартное отклонение, пропорциональное базовому распределению, деленному на квадр.корень из числа выборок.\n", "\n", "Если у вас есть $N$ выборок, среднее значение выборок распределяется по нормальному распределению вокруг истинного среднего со стандартным отклонением $\\sigma/\\sqrt{N}$. \n", "\n", "*Другое* следствие - если перейти от $N_1$ точек данных к $N_2$ точкам данных, можно определить среднее значение $\\sqrt{N_2/N_1}$ более точно. \n", "Увеличение количества данных в 4 раза не дает в 4 раза большей точности, но только удваивает точность." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# получим гистограмму для 1000 значений\n", "data0=[my_get_data(1000).mean() for i in range(1000)]\n", "plt.hist(data0, bins=20)\n", "# настроим ось Х в диапазон 1,5 -- 3,0\n", "plt.xlim(1.5,3)\n", "ti1000=sst.t.interval(0.95, loc=np.mean(data0), df=len(data0)-1, scale=sst.sem(data0)) \n", "## нарисуем верт.линию по значению среднего\n", "plt.axvline(np.mean(data0), label=\"mean\",c='r', ls=\"--\")\n", "plt.axvline(ti1000[0], label=\"mean-\",c='y', ls=\"--\", alpha=0.5)\n", "plt.axvline(ti1000[1], label=\"mean+\",c='y', ls=\"--\", alpha=0.5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1. Распределение выборочных средних приближается к нормальному\n", "2. Ширина (дисперсия) определяется количеством точек, используемых для вычисления каждого выборочного среднего" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n = 1000\n", "data1 = my_get_data(n)\n", "sample_mean = np.mean(data1)\n", "uncert_mean = np.std(data1) / np.sqrt(n)\n", "plt.hist(data1, bins=20)\n", "# получим 95%-доверительный интервал для среднего\n", "ti3=sst.t.interval(0.95, loc=np.mean(data1), df=len(data1)-1, scale=sst.sem(data1))\n", "\n", "## нарисуем верт.линии по значению среднего и доверит.интервала\n", "plt.axvline(np.mean(data1), label=\"mean\",c='r', ls=\"--\")\n", "plt.axvline(ti3[0], label=\"mean-\",c='y', ls=\"--\", alpha=0.5)\n", "plt.axvline(ti3[1], label=\"mean+\",c='y', ls=\"--\", alpha=0.5)\n", "print(f\"Среднее значение совокупности ={sample_mean:.3f} +/- {uncert_mean:.3f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Наличие знака $\\pm$ отражает неопределенность по нормальному распределению." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# VI. Непараметрические критерии\n", "**Что делать, если аналитическое распределение не известно?**\n", "\n", "Рассмотрим пример теста Колмогорова-Смирнова. Обратите внимание, что он предназначен для *непрерывных* распределений, мы используем его для дискретных для того, чтобы понимать пример, но надо следить за состоянием данных, чтобы быть уверенными в результатах, поэтому рассмотрим достаточно \"непрерывное\" дискретное распределение из файла 'loaded__100.txt'." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "r_obs = np.loadtxt(\"loaded__100.txt\")\n", "# np.unique - найдем уникальные значения в массиве\n", "# возвращает отсортированные элементы массива\n", "# Если return_counts =True, также вернет массив счетчиков для каждого элемента, т.е частоты.\n", "\n", "unique1, count_obs = np.unique(r_obs, return_counts=True)\n", "count_obs = count_obs / count_obs.sum() # нормализация\n", "\n", "# для сравнения получим массив равномерно распределенных целых чисел\n", "r_fair = np.random.randint(low=1, high=7, size=(50000, 10)).sum(axis=1)\n", "\n", "# определим набор уникальных значений и их частоту встречаемости\n", "unique2, count_fair = np.unique(r_fair, return_counts=True)\n", "count_fair = count_fair / count_fair.sum() # нормализация\n", "\n", "plt.plot(unique1, count_obs, label=\"мои данные\")\n", "plt.plot(unique2, count_fair, label=\"теоретич.\")\n", "plt.legend(loc=1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# удобно построить вид распределения в коробчатой диаграмме (ящик-с-усами)\n", "plt.boxplot(count_obs, vert=False, patch_artist=True)\n", "print()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Набор данных выглядит довольно беспорядочно.\n", "Если сравнивать по PDF (диаграмме плотности), то приходится объединять результаты в подинтервалы, и трудно получить хорошую диаграмму, если нет большого количества данных.\n", "В таком случае можно построить график распределения (CDF), для которого не нужно формировать данные в гистограмму." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cdf_obs = count_obs.cumsum()\n", "cdf_fair = count_fair.cumsum()\n", "'''\n", "cumsum() — функция в NumPy, которая вычисляет кумулятивные суммы элементов массива по всему массиву или по заданной оси, если указана ось.\n", "Синтаксис функции: numpy.cumsum(array, axis=None, dtype=None, out=None)\n", "Аргументы функции: array — входной массив; axis — ось, по которой вычисляется кумулятивная сумма, если не указано, массив сглаживается.\n", "sa = np.array([[100, 200, 150, 220], [300, 200, 150, 100]])\n", "cumsum = np.cumsum(sa)\n", "# cumsum => ([ 100, 300, 450, 670, 970, 1170, 1320, 1420])\n", "'''\n", "plt.plot(unique1, cdf_obs, label=\"мои данные\")\n", "plt.plot(unique2, cdf_fair, label=\"теоретич.\")\n", "plt.legend(loc=2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Критерий согласия Колмогорова-Смирнова\n", "\n", "*Вызов: kstest (rvs, cdf, args=())*\n", "\n", "где \n", "\n", " rvs --- выборка\n", "\n", " cdf --- функция распределения (сама функция или ее название)\n", "\n", " args --- параметры распределения\n", "\n", "Воспользуемся тестом Колмогорова-Смирнова из scipy.stats для двух выборок." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# отобразим вид наших данных\n", "plt.hist(r_obs, bins=12, density=True, label='histo', color='orange')\n", "# проверка распределений двух выборок\n", "statistic, pvalue = sst.ks_2samp(r_obs, r_fair)\n", "print(f\"Для двух выборок KS-статистика ={statistic:.3f}, p-value ={pvalue:.3f}\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# проверка распределений выборки и равномерного распределения\n", "stats, pval = sst.kstest(r_obs, cdf=sst.uniform.cdf)\n", "print(f\"Для равномерного распределения KS-статистика ={stats:.3f}, p-value ={pval:.3f}\")\n", "\n", "# проверка распределений выборки и нормального распределения\n", "print(sst.kstest(r_obs, sst.norm.cdf))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Подготовим удобную функцию для использования теста Колмогорова-Смирнова для проверки нормальности с визуализацией гистограммы" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def apply_kstest(sample, cdf=sst.norm.cdf):\n", " print(sst.kstest(sample, cdf=cdf))\n", " # ECDF(sample) -- эмпирическая CDF\n", " xx = np.sort(sample)\n", " yy = np.arange(len(xx))/float(len(xx))\n", " grid = np.linspace(-3, 3, 500)\n", "\n", " #plt.figure(figsize=(8, 3.5))\n", " plt.plot(xx, yy, c='red', label='эмпирич.cdf', lw=3)\n", " plt.plot(grid, cdf(grid), c='blue', label='теоретич.cdf', lw=2)\n", " plt.legend()\n", "\n", "apply_kstest(sst.norm.rvs(size=100))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Обратите внимание - Критерий не чувствителен к хвостам распределений. Например, проверим распределение Коши" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "apply_kstest(sst.cauchy.rvs(size=100))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Тест Пирсона $\\chi^2$**\n", "\n", "Критерий согласия Пирсона (хи-квадрат)\n", "\n", "𝖧0 : Выборка из некоторого класса распределений\n", "\n", "𝖧1 : 𝖧0 не верна\n", "\n", "$$ \\chi^2 = \\sum_{i=1}^{n} \\frac{(C_i - E_i)^2}{E_i} $$\n", "\n", "Что значит статистика $\\chi^2$ ? Можно сконвертировать результат в вероятность по $\\chi^2$ распределению" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# chi2 - Непрерывная случайная величина хи-квадрат\n", "df = 55 # примем количество степеней свободы =55\n", "# Вычислим первые четыре момента выборки \n", "mean, var, skew, kurt = sst.chi2.stats(df, moments='mvsk')\n", "\n", "x = np.linspace(sst.chi2.ppf(0.01, df),\n", " sst.chi2.ppf(0.99, df),\n", " 100)\n", "# Отобразим функцию плотности вероятности (pdf -Probability Density Function)\n", "# pdf(x) — значение плотности в точке x \n", "plt.plot(x, sst.chi2.pdf(x, df), 'r-', lw=5, alpha=0.6, label='chi2 pdf')\n", "rv = sst.chi2(df)\n", "\n", "# Проверим точность cdf и ppf\n", "plt.plot(x, rv.pdf(x), 'k-', lw=2, label='chi2 pdf')\n", "vals = sst.chi2.ppf([0.001, 0.5, 0.999], df)\n", "np.allclose([0.001, 0.5, 0.999], sst.chi2.cdf(vals, df))\n", "\n", "# Генерация случайных чисел\n", "r = sst.chi2.rvs(df, size=1000)\n", "\n", "# сравним гистограмму\n", "plt.hist(r, density=True, bins='auto', histtype='stepfilled', alpha=0.2)\n", "plt.xlim([x[0], x[-1]])\n", "plt.legend(loc='best', frameon=False)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "## используем функцию chi2 из фреймворка scipy.stats\n", "chi2s = np.linspace(0, 15, 500)\n", "prob = sst.chi2.pdf(chi2s, 5)\n", "\n", "chi2_val = 7 # число выбрано просто для примера!\n", "## нарисуем график распределения\n", "plt.plot(chi2s, prob, label=\"распределение\")\n", "## нарисуем верт.линию по значению статистики\n", "plt.axvline(chi2_val, label=\"$\\\\chi2$\", ls=\"--\")\n", "## заштрихуем критич.область\n", "plt.fill_between(chi2s, prob, 0, where=(chi2s>=chi2_val), alpha=0.1)\n", "## выведем на диаграмму легенду\n", "plt.legend()\n", "# sf - Функция выживания, она определяется как (1 - cdf) // Survival function\n", "## выведем на консоль\n", "print(f\"Значение p-value = {sst.chi2.sf(chi2_val, 5):.3f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Вызов: chisquare (f_obs, f_exp=None, ddof=0)*\n", "\n", "где\n", "\n", " f_obs --- число элементов выборки, попавших в каждый из интервалов\n", "\n", " f_exp --- ожидаемое число (по умолчанию равномерное)\n", "\n", " ddof --- поправка на число степеней свободы. Статистика асимптотически будет иметь распределение 𝑘 − 1 − 𝑑𝑑𝑜𝑓, где 𝑘 -- число интервалов. \n", "\n", "Несколько тестов могут быть действительны. Будьте осторожны при поиске нужного способа проверки гипотезы - это называется \"поиск значимости\". Чем конкретнее ваш \"вопрос\", тем более мощный тест вы можете использовать.\n", "\n", "В этом примере используется одностороннее распределение, которое проверяет вопрос, не является ли наше наблюдаемое распределение *слишком отличным* от базового.\n", "Можно также использовать двустороннее распределение, которое проверяет, не является ли наше наблюдаемое распределение *слишком похожим* на базовое.\n", "\n", "Для получения своих данных выгрузите из нашей модели Excel семинара 2 столбец W в формате txt (файл- 12345.prn)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from collections import Counter\n", "# загрузим для проверки свои данные- в формате строк\n", "s_obs = np.loadtxt(\"12345.prn\",dtype=\"str\")\n", "# теперь строки преобразуем в массив numpy, конвертируя числа с десятичной запятой вместо десятичной точки (при необходимости!)\n", "s_ob = np.array([float(so.replace(',', '.')) for so in s_obs])\n", "# выведем полученный результат с основными параметрами\n", "print(f\"прочитано {len(s_obs)} строк |среднее значение = {s_ob.mean():.3f} |стандарт.отклонение = {s_ob.std():.3f} |максим.значение = {s_ob.max():.3f}\")\n", "\n", "desc_ob = sst.describe(s_ob, axis=0, ddof=1, bias=True, nan_policy='propagate')\n", "print(desc_ob)\n", "\n", "# stats(df, loc=0, scale=1, moments=’mv’)\n", "# Mean(‘m’), variance(‘v’), skew(‘s’), kurtosis(‘k’)\n", "f_obs = np.array(list(Counter(s_ob).values()))\n", "s_hist = np.histogram(s_ob, bins=10)\n", "s_fair = np.random.normal(s_ob.mean(), s_ob.std(), size=(100, 1))\n", "f_hist = np.histogram(s_fair, bins=10)\n", "#\n", "plt.hist(s_fair, bins=10)\n", "plt.hist(s_ob, bins=10)\n", "#\n", "chisq1, pv1 = sst.chisquare(s_hist[0], f_hist[0])\n", "print(f\"Значение статистики chi2 = {chisq1:.3f} и значение p-value = {pv1:.3f}\")\n", "if pv1 < 0.05:\n", " print(\"-Гипотезу H0 можно отклонить\")\n", "else:\n", " print(\"-Гипотезу H0 не отклоняем\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Примените тест Колмогорова к своим данным (возможно, немного смасштабированным для лучшей визуализации)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "## тест Колмогорова-Смирнова распределения своей выборки и равномерного распределения" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__График QQ__\n", "\n", "У вероятностного графика в математической статистике данные выводят на график относительно теоретического распределения таким образом, чтобы точки образовывали примерно прямую линию. Отклонения от этой прямой указывают на отклонения от заданного распределения.\n", "\n", "График QQ («Q» обозначает квантиль ) — вероятностный график, представляющий собой графический метод сравнения двух распределений вероятностей путем расставления их квантилей друг против друга. Для этого надо задать набор интервалов для квантилей. Точке (х, у) на графике соответствует один из квантилей второго распределения (у-координата), построенная вместе с аналогичным квантилем первого распределения (х-координата). Таким образом, линия является параметрической кривой с параметром, равным количеству интервалов квантилей.\n", "\n", "Если два сравниваемых распределения схожи, точки в графике QQ будут примерно лежать на прямой у = х. Если распределения линейно связаны, точки в графике QQ будут примерно лежать на одной прямой, но не обязательно на прямой у = х. График QQ также может быть использован в качестве графического средства оценки параметров в масштабе семейства распределений.\n", "\n", "График QQ используется для сравнения формы распределения, обеспечивая графическое представление о свойствах, таких как местоположение, масштаб и асимметрия в двух распределениях. График QQ может быть использован для сравнения с набором данных или с теоретическим распределением. Использование графика QQ для сравнения двух наборов данных можно рассматривать как непараметрический подход к сравнению их распределений. \n", "\n", "График QQ — более мощный подход, чем метод сравнения гистограмм двух образцов, но требует больше умений для интерпретации. График QQ, как правило, используется для сравнения данных с теоретической моделью. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def QQ_plot(data):\n", " # сортировка по возрастанию\n", " yo = np.sort(data)\n", " \n", " # получить среднее и станд.отклонение\n", " mean, std = np.mean(yo), np.std(yo)\n", " \n", " # получить множество квантилей нормал.распределения\n", " # ppf(q) — q-й квантиль (Percent Point Function) - Inverse CDF\n", " ppf = sst.norm(loc=mean, scale=std).ppf\n", " xo = [ppf( i/(len(yo)+2) ) for i in range(1,len(yo)+1)]\n", "\n", " # построить точечный график Q-Q\n", " plt.scatter(xo, yo)\n", " \n", " # построить диагональную линию\n", " diag = np.linspace(np.min([xo,yo]), np.max([xo,yo]), 1000)\n", " plt.plot(diag, diag, c='red', linestyle='--')\n", "\n", " # график Q-Q в квадратной форме\n", " plt.gca().set_aspect('equal')\n", " \n", " # добавить подписи данных\n", " plt.xlabel('квантили нормал.распр.')\n", " plt.ylabel('квантили моих данных')\n", "\n", "# сделаем для теста произвольные нормальные данные - через numpy или scipy\n", "nx = np.random.normal(loc=0.004, scale=0.002, size=100)\n", "#nx = sst.norm(loc=0.004, scale=0.002).rvs(size=100)\n", "QQ_plot(nx)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Постройте QQ-график для своих данных." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# используем тест Q-Q для своих данных" ] } ], "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.12.7" } }, "nbformat": 4, "nbformat_minor": 2 }