ДЗ3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Задание 3.1\n",
    "\n",
    "## Гипергеометрическое распределени\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Свойства выборочного среднего и выборочной дисперсии.\n",
    "\n",
    "Пусть $X_1$, $X_2$, ..., $X_n$ - независимые одинаково распределённые случайные величины. <br>\n",
    "$\\bar{X} = \\frac{1}{n} \\sum_{i=1}^{n}X_i$ - выборочное среднее, $S_n^2 = \\frac{1}{n} \\sum_{i=1}^{n}(X_i - \\bar{X})^2$ - выборочная дисперсия <br>\n",
    "<br>\n",
    "Тогда $E[\\bar{X}] = E[\\frac{1}{n} \\sum_{i=1}^{n}X_i] = \\frac{1}{n} \\sum_{i=1}^{n}E[X_i] = \\frac{1}{n} \\sum_{i=1}^{n}E[X_1] = \\frac{n}{n} E[X_1] = E[X_1]$ <br>\n",
    "Что является несмещённой оценкой матожидания исследуемого распределения. <br><br>\n",
    "Поскольку $S_n^2$ оценивает дисперсию исследуемого распределения, то можно рассмотреть выборочную дисперсию для\n",
    "величины $Y_i= X_i - E[X_1], i=1...n$<br>\n",
    "Тогда $E[Y_i]= E[\\bar{Y}] = 0$, $D[X_i] = D[Y_i] = E[Y_i^2]$. <br>\n",
    "$ S_n^2 = \\frac{1}{n} \\sum_{i=1}^{n}(Y_i - \\bar{Y})^2 =  \\frac{1}{n} \\sum_{i=1}^{n} Y_i^2 - \\frac{2}{n} \\sum_{i=1}^{n} Y_i \\bar{Y} + \\frac{1}{n} \\sum_{i=1}^{n} \\bar{Y}^2 = \\frac{1}{n} \\sum_{i=1}^{n} Y_i^2 - 2 \\bar{Y} \\frac{1}{n}\\sum_{i=1}^{n} Y_i   + \\bar{Y}^2 = \\frac{1}{n} \\sum_{i=1}^{n} Y_i^2 - 2 \\bar{Y}^2 + \\bar{Y}^2 =$ <br> $ = \\frac{1}{n} \\sum_{i=1}^{n} Y_i^2 - \\bar{Y}^2  $ <br>\n",
    "$E[S_n^2] = \\frac{1}{n} \\sum_{i=1}^{n}E[Y_i^2] - E[\\bar{Y}^2] = D[Y_1] - E[\\bar{Y}^2] $ <br>\n",
    "Заметим, что $E[Y_i Y_j] = 0$ для $i\\neq j$, так как случайные величины являются независимыми <br>\n",
    "$ E[\\bar{Y}^2] = \\frac{1}{n^2} E[(\\sum_{i=1}^{n}Y_i)^2] = \\frac{1}{n^2} \\sum_{i=1}^{n} E[Y_i^2] + \\sum_{i\\neq j} E[Y_i Y_j] = \n",
    "\\frac{1}{n} E[Y_1^2] = \\frac{1}{n} D[Y_1] $ <br> \n",
    "Имеем:  $E[S_n^2] =  D[Y_1] - \\frac{1}{n} D[Y_1] = \\frac{n-1}{n} D[Y_1] $ <br><br>\n",
    "Эта оценка является состоятельной, но смещённой.\n",
    "<br><br>\n",
    "Несмещённой состоятельной оценкой выборочной дисперсии будет  $S^2 = \\frac{n}{n-1} S_n^2 = \\frac{1}{n-1} \\sum_{i=1}^{n}(Y_i - \\bar{Y})^2 $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Зададим гипергеометрическое распределение $HG(200, 70, 10)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import hypergeom\n",
    "\n",
    "volumes = [5, 10, 100, 1000, 10**5]\n",
    "count_samples = 5\n",
    "\n",
    "N, M, n = (200, 70, 10)\n",
    "hg = hypergeom(N, n, M)\n",
    "\n",
    "np.random.seed(1000)\n",
    "data = {v: np.array([hg.rvs(size=v) for _ in range(count_samples)]) for v in volumes}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Подсчёт выборочного среднего и выборочной дисперсии\n",
    "\n",
    "Отметим, что для $\\xi \\sim HG (200, 70, 10)$: <br>\n",
    "$E[\\xi] = 3.5  $ <br>\n",
    "$D[\\xi] = \\frac{1729}{796} \\approx 2.17211  $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     n = 5, Выборка №1, выборочное среднее = 4.0, выборочная дисперсия = 1.2, несмещённая выборочная дисперсия = 1.5\n",
      "     n = 5, Выборка №2, выборочное среднее = 2.8, выборочная дисперсия = 0.56, несмещённая выборочная дисперсия = 0.7\n",
      "     n = 5, Выборка №3, выборочное среднее = 3.2, выборочная дисперсия = 1.36, несмещённая выборочная дисперсия = 1.7\n",
      "     n = 5, Выборка №4, выборочное среднее = 4.2, выборочная дисперсия = 2.16, несмещённая выборочная дисперсия = 2.7\n",
      "     n = 5, Выборка №5, выборочное среднее = 3.2, выборочная дисперсия = 0.56, несмещённая выборочная дисперсия = 0.7\n",
      "\n",
      "    n = 10, Выборка №1, выборочное среднее = 4.1, выборочная дисперсия = 2.09, несмещённая выборочная дисперсия = 2.322\n",
      "    n = 10, Выборка №2, выборочное среднее = 2.3, выборочная дисперсия = 2.21, несмещённая выборочная дисперсия = 2.456\n",
      "    n = 10, Выборка №3, выборочное среднее = 2.8, выборочная дисперсия = 2.16, несмещённая выборочная дисперсия = 2.4\n",
      "    n = 10, Выборка №4, выборочное среднее = 3.6, выборочная дисперсия = 2.44, несмещённая выборочная дисперсия = 2.711\n",
      "    n = 10, Выборка №5, выборочное среднее = 3.4, выборочная дисперсия = 2.84, несмещённая выборочная дисперсия = 3.156\n",
      "\n",
      "   n = 100, Выборка №1, выборочное среднее = 3.59, выборочная дисперсия = 2.02, несмещённая выборочная дисперсия = 2.042\n",
      "   n = 100, Выборка №2, выборочное среднее = 3.3, выборочная дисперсия = 1.79, несмещённая выборочная дисперсия = 1.808\n",
      "   n = 100, Выборка №3, выборочное среднее = 3.82, выборочная дисперсия = 2.41, несмещённая выборочная дисперсия = 2.432\n",
      "   n = 100, Выборка №4, выборочное среднее = 3.44, выборочная дисперсия = 2.17, несмещённая выборочная дисперсия = 2.188\n",
      "   n = 100, Выборка №5, выборочное среднее = 3.77, выборочная дисперсия = 2.2, несмещённая выборочная дисперсия = 2.219\n",
      "\n",
      "  n = 1000, Выборка №1, выборочное среднее = 3.51, выборочная дисперсия = 2.15, несмещённая выборочная дисперсия = 2.148\n",
      "  n = 1000, Выборка №2, выборочное среднее = 3.47, выборочная дисперсия = 2.13, несмещённая выборочная дисперсия = 2.133\n",
      "  n = 1000, Выборка №3, выборочное среднее = 3.41, выборочная дисперсия = 2.13, несмещённая выборочная дисперсия = 2.132\n",
      "  n = 1000, Выборка №4, выборочное среднее = 3.53, выборочная дисперсия = 2.31, несмещённая выборочная дисперсия = 2.317\n",
      "  n = 1000, Выборка №5, выборочное среднее = 3.51, выборочная дисперсия = 2.26, несмещённая выборочная дисперсия = 2.266\n",
      "\n",
      "n = 100000, Выборка №1, выборочное среднее = 3.5, выборочная дисперсия = 2.17, несмещённая выборочная дисперсия = 2.172\n",
      "n = 100000, Выборка №2, выборочное среднее = 3.5, выборочная дисперсия = 2.16, несмещённая выборочная дисперсия = 2.158\n",
      "n = 100000, Выборка №3, выборочное среднее = 3.49, выборочная дисперсия = 2.17, несмещённая выборочная дисперсия = 2.168\n",
      "n = 100000, Выборка №4, выборочное среднее = 3.5, выборочная дисперсия = 2.16, несмещённая выборочная дисперсия = 2.156\n",
      "n = 100000, Выборка №5, выборочное среднее = 3.49, выборочная дисперсия = 2.19, несмещённая выборочная дисперсия = 2.187\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for volume in volumes:\n",
    "    for choice in range(count_samples):\n",
    "        print(f\"{'n = '+str(volume):>10},\",\n",
    "              f\"Выборка №{choice+1}, выборочное среднее = {np.mean(data[volume][choice]):.3},\",\n",
    "              f\"выборочная дисперсия = {np.var(data[volume][choice]):.3},\",\n",
    "              f\"несмещённая выборочная дисперсия = {np.var(data[volume][choice], ddof=1):.4}\")\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как видно из при увеличении n значения оценок сходятся к истинным значениям параметров"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Гамма распределение\n",
    "\n",
    "### Зададим гамма распределение $Г(2, 1/2)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import gamma\n",
    "\n",
    "volumes = [5, 10, 100, 1000, 10**5]\n",
    "count_samples = 5\n",
    "\n",
    "a, th = 2, 1/2 \n",
    "g = gamma(a, scale=1/th)\n",
    "\n",
    "np.random.seed(1000)\n",
    "data = {v: np.array([g.rvs(size=v) for _ in range(count_samples)]) for v in volumes}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Подсчёт выборочного среднего и выборочной дисперсии\n",
    "\n",
    "Отметим, что для $\\xi \\sim Г(2, 1/2)$: <br>\n",
    "$E[\\xi] = 4  $ <br>\n",
    "$D[\\xi] = 8  $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     n = 5, Выборка №1, выборочное среднее = 3.61, выборочная дисперсия = 1.63, несмещённая выборочная дисперсия = 2.039\n",
      "     n = 5, Выборка №2, выборочное среднее = 2.85, выборочная дисперсия = 2.14, несмещённая выборочная дисперсия = 2.672\n",
      "     n = 5, Выборка №3, выборочное среднее = 4.42, выборочная дисперсия = 5.94, несмещённая выборочная дисперсия = 7.429\n",
      "     n = 5, Выборка №4, выборочное среднее = 4.59, выборочная дисперсия = 11.0, несмещённая выборочная дисперсия = 13.75\n",
      "     n = 5, Выборка №5, выборочное среднее = 2.6, выборочная дисперсия = 2.3, несмещённая выборочная дисперсия = 2.876\n",
      "\n",
      "    n = 10, Выборка №1, выборочное среднее = 1.57, выборочная дисперсия = 1.3, несмещённая выборочная дисперсия = 1.445\n",
      "    n = 10, Выборка №2, выборочное среднее = 4.23, выборочная дисперсия = 12.4, несмещённая выборочная дисперсия = 13.82\n",
      "    n = 10, Выборка №3, выборочное среднее = 4.65, выборочная дисперсия = 7.92, несмещённая выборочная дисперсия = 8.796\n",
      "    n = 10, Выборка №4, выборочное среднее = 3.32, выборочная дисперсия = 4.33, несмещённая выборочная дисперсия = 4.81\n",
      "    n = 10, Выборка №5, выборочное среднее = 2.61, выборочная дисперсия = 1.14, несмещённая выборочная дисперсия = 1.265\n",
      "\n",
      "   n = 100, Выборка №1, выборочное среднее = 3.32, выборочная дисперсия = 5.64, несмещённая выборочная дисперсия = 5.694\n",
      "   n = 100, Выборка №2, выборочное среднее = 4.58, выборочная дисперсия = 7.22, несмещённая выборочная дисперсия = 7.298\n",
      "   n = 100, Выборка №3, выборочное среднее = 4.43, выборочная дисперсия = 9.29, несмещённая выборочная дисперсия = 9.384\n",
      "   n = 100, Выборка №4, выборочное среднее = 3.63, выборочная дисперсия = 6.48, несмещённая выборочная дисперсия = 6.542\n",
      "   n = 100, Выборка №5, выборочное среднее = 3.77, выборочная дисперсия = 6.27, несмещённая выборочная дисперсия = 6.338\n",
      "\n",
      "  n = 1000, Выборка №1, выборочное среднее = 3.91, выборочная дисперсия = 8.21, несмещённая выборочная дисперсия = 8.22\n",
      "  n = 1000, Выборка №2, выборочное среднее = 4.0, выборочная дисперсия = 7.97, несмещённая выборочная дисперсия = 7.979\n",
      "  n = 1000, Выборка №3, выборочное среднее = 4.06, выборочная дисперсия = 8.67, несмещённая выборочная дисперсия = 8.681\n",
      "  n = 1000, Выборка №4, выборочное среднее = 4.08, выборочная дисперсия = 8.66, несмещённая выборочная дисперсия = 8.673\n",
      "  n = 1000, Выборка №5, выборочное среднее = 3.91, выборочная дисперсия = 6.97, несмещённая выборочная дисперсия = 6.979\n",
      "\n",
      "n = 100000, Выборка №1, выборочное среднее = 4.0, выборочная дисперсия = 8.03, несмещённая выборочная дисперсия = 8.027\n",
      "n = 100000, Выборка №2, выборочное среднее = 3.99, выборочная дисперсия = 7.91, несмещённая выборочная дисперсия = 7.906\n",
      "n = 100000, Выборка №3, выборочное среднее = 4.01, выборочная дисперсия = 7.95, несмещённая выборочная дисперсия = 7.948\n",
      "n = 100000, Выборка №4, выборочное среднее = 4.03, выборочная дисперсия = 8.15, несмещённая выборочная дисперсия = 8.151\n",
      "n = 100000, Выборка №5, выборочное среднее = 3.99, выборочная дисперсия = 7.95, несмещённая выборочная дисперсия = 7.951\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for volume in volumes:\n",
    "    for choice in range(count_samples):\n",
    "        print(f\"{'n = '+str(volume):>10},\",\n",
    "              f\"Выборка №{choice+1}, выборочное среднее = {np.mean(data[volume][choice]):.3},\",\n",
    "              f\"выборочная дисперсия = {np.var(data[volume][choice]):.3},\",\n",
    "              f\"несмещённая выборочная дисперсия = {np.var(data[volume][choice], ddof=1):.4}\")\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как видно из при увеличении n значения оценок сходятся к истинным значениям параметров"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Задание 3.2\n",
    "\n",
    "## Гипергеометрическое распределение\n",
    "\n",
    "### Теоретическое построение оценок"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Рассмотрим задачу оценки численности замкнутой популяции животных.\n",
    "Тогда будем считать, что в какой-то момент времени  $M$ животных были помечены бирками, и в настоящий момент имеется $m$ выборок по $n$ наблюдейний в каждом.<br>\n",
    "$x_1, x_2, ... , x_m \\sim HG(\\theta, M, n)$, где $\\theta$ - неизвестная численность популяции. <br>\n",
    "Найдём оценку методом моментов. <br>\n",
    "Статистика $T(\\vec{X}) = \\bar{x}$ является смещённой оценокой для $\\theta$, так как $E[T]=E[x_1]=\\frac{M n}{\\theta}$ <br>\n",
    "По методу моментов оценкой неизвесного параметра будет $\\hat{\\theta}=\\frac{M n}{\\bar{x}}$. Эта оценка является состоятельной. <br><br>\n",
    "\n",
    "Теперь найдём оценку для численности популяции методом максимального правдоподобия. <br>\n",
    "Для этого предположим, что $\\frac{M}{\\theta} << 1$ и $n >> 1$, тогда гиепергеометрическое распределение приближается распределением Пуассона $\\Pi(\\frac{nM}{\\theta}) $. [1]<br>\n",
    "Логарифмическая функция правдоподобия будет иметь вид: <br>\n",
    "$ln(L(\\vec{X}, \\theta)) = \\sum_{i=1}^{m}ln(\\frac{exp\\{-\\frac{nM}{\\theta}\\}(\\frac{nM}{\\theta})^{x_i}}{x_i!})= -m\\frac{nM}{\\theta}-ln(\\theta)\\sum_{i=1}^{m}x_i  + ln(nM)\\sum_{i=1}^{m}x_i - \\sum_{i=1}^{m}ln(x_i!)$<br><br>\n",
    "Найдём максимум этой функции:<br>\n",
    "$$\\frac{\\partial ln L}{\\partial \\theta} = 0$$\n",
    "$$m \\frac{nM}{\\hat\\theta^2}-\\frac{1}{\\hat\\theta}\\sum_{i=1}^{m}x_i=0 \\space \\space \\Leftrightarrow \\space\\space \n",
    "nM = \\hat\\theta\\frac{1}{m}\\sum_{i=1}^{m}x_i \\space \\space \\Leftrightarrow \\space\\space \\hat{\\theta}=\\frac{nM}{\\bar{x}}$$<br>\n",
    "Оценка методом максимального правдоподобия точно такой же, как и оценка методом моментов. <br>\n",
    "\n",
    "<br>\n",
    "Стоит отметить, что исходная модель не является регулярной, однако после приближения распределением Пуассона, регулярность появилась. <br>\n",
    "Продолжим исследавать модель и найдём эффективную оценку и параметрическую функцию, которую она оценивает. <br>\n",
    "$$V = \\frac{\\partial ln L}{\\partial \\theta} = m \\frac{nM}{\\theta^2}-\\frac{1}{\\theta}\\sum_{i=1}^{m}x_i \\space \\space \\Leftrightarrow \\space\\space \\frac{\\theta}{m nM} V = \\frac{1}{\\theta} - \\frac{\\bar{x}}{nM} $$<br>\n",
    "Таким образом, $T(\\vec{X}) = \\frac{\\bar{x}}{nM}$ является эффективной оценкой для $\\frac{1}{\\theta}$ <br>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[1] https://stats.stackexchange.com/questions/73159/verification-of-poisson-approximation-to-hypergeometric-distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Проверка оценок на данных"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     n = 5, Выборка №1, оцениваемый параметр = 175,\n",
      "     n = 5, Выборка №2, оцениваемый параметр = 250,\n",
      "     n = 5, Выборка №3, оцениваемый параметр = 218,\n",
      "     n = 5, Выборка №4, оцениваемый параметр = 166,\n",
      "     n = 5, Выборка №5, оцениваемый параметр = 218,\n",
      "\n",
      "    n = 10, Выборка №1, оцениваемый параметр = 170,\n",
      "    n = 10, Выборка №2, оцениваемый параметр = 304,\n",
      "    n = 10, Выборка №3, оцениваемый параметр = 250,\n",
      "    n = 10, Выборка №4, оцениваемый параметр = 194,\n",
      "    n = 10, Выборка №5, оцениваемый параметр = 205,\n",
      "\n",
      "   n = 100, Выборка №1, оцениваемый параметр = 194,\n",
      "   n = 100, Выборка №2, оцениваемый параметр = 212,\n",
      "   n = 100, Выборка №3, оцениваемый параметр = 183,\n",
      "   n = 100, Выборка №4, оцениваемый параметр = 203,\n",
      "   n = 100, Выборка №5, оцениваемый параметр = 185,\n",
      "\n",
      "  n = 1000, Выборка №1, оцениваемый параметр = 199,\n",
      "  n = 1000, Выборка №2, оцениваемый параметр = 201,\n",
      "  n = 1000, Выборка №3, оцениваемый параметр = 205,\n",
      "  n = 1000, Выборка №4, оцениваемый параметр = 198,\n",
      "  n = 1000, Выборка №5, оцениваемый параметр = 199,\n",
      "\n",
      "n = 100000, Выборка №1, оцениваемый параметр = 200,\n",
      "n = 100000, Выборка №2, оцениваемый параметр = 199,\n",
      "n = 100000, Выборка №3, оцениваемый параметр = 200,\n",
      "n = 100000, Выборка №4, оцениваемый параметр = 200,\n",
      "n = 100000, Выборка №5, оцениваемый параметр = 200,\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for volume in volumes:\n",
    "    for choice in range(count_samples):\n",
    "        print(f\"{'n = '+str(volume):>10},\",\n",
    "              f\"Выборка №{choice+1}, оцениваемый параметр = {int(n*M/np.mean(data[volume][choice])):},\")\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Гамма распределение\n",
    "\n",
    "### Теоретическое построение оценок"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В реальной ситуации интерес представляет выяспенение обоих параметров гамма распределения. Сформулируем задачу. <br>\n",
    "$x_1, x_2, ... , x_n \\sim \\Gamma(\\theta_1, \\theta_2)$, где $\\vec{\\theta} = (\\theta_1, \\theta_2)$ - неизвестные параметры. <br>\n",
    "Напишем систему уравнений по методу момнтов.<br>\n",
    "$$\\begin{cases} \\frac{\\hat{\\theta_1}}{\\hat{\\theta_2}} = \\bar{x} \\\\ \\frac{\\hat{\\theta_1}}{\\hat{\\theta_2}^2} = s^2 \\end{cases} \\space \\space \\Leftrightarrow \\space\\space  \\begin{cases}  \\hat{\\theta_1}= \\frac{\\bar{x}^2}{s^2}  \\\\ \\hat{\\theta_2} =  \\frac{\\bar{x}}{s^2} \\end{cases},  $$\n",
    "где $\\bar{x}$ - выборочное среднее, $s^2$ - несмещённая выборочная дисперсия. <br>\n",
    "Полученные оценки являются состоятельными. Вопрос об их смещённости остаётся открытым. <br>\n",
    "Проверим эту оценку на данных."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     n = 5, Выборка №1, первый параметр = 6.4,второй параметр = 1.772\n",
      "     n = 5, Выборка №2, первый параметр = 3.04,второй параметр = 1.067\n",
      "     n = 5, Выборка №3, первый параметр = 2.63,второй параметр = 0.5953\n",
      "     n = 5, Выборка №4, первый параметр = 1.53,второй параметр = 0.3335\n",
      "     n = 5, Выборка №5, первый параметр = 2.36,второй параметр = 0.9055\n",
      "\n",
      "    n = 10, Выборка №1, первый параметр = 1.71,второй параметр = 1.089\n",
      "    n = 10, Выборка №2, первый параметр = 1.3,второй параметр = 0.3063\n",
      "    n = 10, Выборка №3, первый параметр = 2.45,второй параметр = 0.5282\n",
      "    n = 10, Выборка №4, первый параметр = 2.29,второй параметр = 0.6905\n",
      "    n = 10, Выборка №5, первый параметр = 5.39,второй параметр = 2.065\n",
      "\n",
      "   n = 100, Выборка №1, первый параметр = 1.93,второй параметр = 0.5827\n",
      "   n = 100, Выборка №2, первый параметр = 2.87,второй параметр = 0.627\n",
      "   n = 100, Выборка №3, первый параметр = 2.09,второй параметр = 0.4722\n",
      "   n = 100, Выборка №4, первый параметр = 2.02,второй параметр = 0.5556\n",
      "   n = 100, Выборка №5, первый параметр = 2.24,второй параметр = 0.5946\n",
      "\n",
      "  n = 1000, Выборка №1, первый параметр = 1.86,второй параметр = 0.4753\n",
      "  n = 1000, Выборка №2, первый параметр = 2.01,второй параметр = 0.5016\n",
      "  n = 1000, Выборка №3, первый параметр = 1.9,второй параметр = 0.4675\n",
      "  n = 1000, Выборка №4, первый параметр = 1.92,второй параметр = 0.4701\n",
      "  n = 1000, Выборка №5, первый параметр = 2.19,второй параметр = 0.5602\n",
      "\n",
      "n = 100000, Выборка №1, первый параметр = 1.99,второй параметр = 0.498\n",
      "n = 100000, Выборка №2, первый параметр = 2.01,второй параметр = 0.5048\n",
      "n = 100000, Выборка №3, первый параметр = 2.03,второй параметр = 0.5048\n",
      "n = 100000, Выборка №4, первый параметр = 1.99,второй параметр = 0.4942\n",
      "n = 100000, Выборка №5, первый параметр = 2.0,второй параметр = 0.5017\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for volume in volumes:\n",
    "    for choice in range(count_samples):\n",
    "        x_sr = np.mean(data[volume][choice])\n",
    "        x_d = np.var(data[volume][choice], ddof=1)\n",
    "        print(f\"{'n = '+str(volume):>10},\",\n",
    "              f\"Выборка №{choice+1}, первый параметр = {x_sr**2 / x_d:.3},\"\n",
    "              f\"второй параметр = {x_sr / x_d:.4}\")\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видна явная состоятельность оценки"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Теперь попробуем найти оценки для гамма распределения с 2-мя неизвестными параметрами методом максимального правдоподобия.<br>\n",
    "Логарифмическая функция правдоподобия будет иметь вид: <br>\n",
    "$ln(L(\\vec{X}, \\vec{\\theta})) = \\sum_{i=1}^{n}((\\theta_1-1)ln(x_i) + \\theta_1 ln(\\theta_2) - \\theta_2 x_i - ln(\\Gamma(\\theta_1))) = (\\theta_1-1)\\sum_{i=1}^{n}ln(x_i) + n \\theta_1 ln(\\theta_2) - \\theta_2 \\sum_{i=1}^{n}x_i - n\\space \n",
    "ln(\\Gamma(\\theta_1))$<br>\n",
    "Отметим несколько фактов. Заметим, что модель является регулярной и экспоненциальной. <br>\n",
    "Достаточными статистиками: <br>\n",
    "<ul>\n",
    " <li>Для $\\theta_1$ ялвяется $T_1(\\vec{X}) = \\sum_{i=1}^{n}ln(x_i)$ </li>\n",
    " <li>Для $\\theta_2$ является $T_1(\\vec{X}) = \\sum_{i=1}^{n}x_i$ </li>\n",
    "</ul>\n",
    "Составим систему урванений для метода максимального правдоподобия. <br>\n",
    "$$\\begin{cases} \\frac{\\partial ln L}{\\partial \\theta_1} = 0 \\\\ \\frac{\\partial ln L}{\\partial \\theta_2} = 0 \\end{cases} \\space \\space \\Leftrightarrow \\space\\space  \\begin{cases} \\sum_{i=1}^{n}ln(x_i) + n \\space ln(\\theta_2) - n \\space \\psi(\\theta_1) = 0  \\\\ n \\theta_1 \\frac{1}{\\theta_2} - \\sum_{i=1}^{n}x_i = 0 \\end{cases} \\space \\space \\Leftrightarrow \\space\\space  \\begin{cases}  \\theta_2 =  \\frac{\\theta_1}{\\bar{X}}\\\\ \\overline{ln(X)} + ln(\\theta_1) - ln(\\bar{X}) -  \\psi(\\theta_1) = 0 \\end{cases}  , $$\n",
    "где $\\psi(\\theta_1) = (ln(\\Gamma(\\theta_1))'$ - диагамма, $\\bar{X} = \\frac{1}{n}\\sum_{i=1}^{n}x_i$, $\\overline{ln(X)} = \\frac{1}{n}\\sum_{i=1}^{n}ln(x_i)$ <br>\n",
    "Второе уравнение является трансцендентным уравнением и поэтому можно решить только численными методами. <br>\n",
    "Проверим оценку на данных."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     n = 5, Выборка №1, первый параметр = 6.58,второй параметр = 1.82\n",
      "     n = 5, Выборка №2, первый параметр = 4.0,второй параметр = 1.4\n",
      "     n = 5, Выборка №3, первый параметр = 2.75,второй параметр = 0.621\n",
      "     n = 5, Выборка №4, первый параметр = 2.83,второй параметр = 0.618\n",
      "     n = 5, Выборка №5, первый параметр = 3.33,второй параметр = 1.28\n",
      "\n",
      "    n = 10, Выборка №1, первый параметр = 1.92,второй параметр = 1.22\n",
      "    n = 10, Выборка №2, первый параметр = 1.45,второй параметр = 0.343\n",
      "    n = 10, Выборка №3, первый параметр = 2.49,второй параметр = 0.536\n",
      "    n = 10, Выборка №4, первый параметр = 1.82,второй параметр = 0.547\n",
      "    n = 10, Выборка №5, первый параметр = 5.23,второй параметр = 2.0\n",
      "\n",
      "   n = 100, Выборка №1, первый параметр = 2.0,второй параметр = 0.603\n",
      "   n = 100, Выборка №2, первый параметр = 2.54,второй параметр = 0.555\n",
      "   n = 100, Выборка №3, первый параметр = 1.85,второй параметр = 0.417\n",
      "   n = 100, Выборка №4, первый параметр = 2.26,второй параметр = 0.622\n",
      "   n = 100, Выборка №5, первый параметр = 2.12,второй параметр = 0.562\n",
      "\n",
      "  n = 1000, Выборка №1, первый параметр = 1.85,второй параметр = 0.474\n",
      "  n = 1000, Выборка №2, первый параметр = 1.99,второй параметр = 0.498\n",
      "  n = 1000, Выборка №3, первый параметр = 1.94,второй параметр = 0.478\n",
      "  n = 1000, Выборка №4, первый параметр = 1.91,второй параметр = 0.468\n",
      "  n = 1000, Выборка №5, первый параметр = 2.07,второй параметр = 0.528\n",
      "\n",
      "n = 100000, Выборка №1, первый параметр = 1.99,второй параметр = 0.499\n",
      "n = 100000, Выборка №2, первый параметр = 2.0,второй параметр = 0.5\n",
      "n = 100000, Выборка №3, первый параметр = 2.01,второй параметр = 0.5\n",
      "n = 100000, Выборка №4, первый параметр = 1.99,второй параметр = 0.494\n",
      "n = 100000, Выборка №5, первый параметр = 2.01,второй параметр = 0.503\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from scipy.optimize import bisect\n",
    "from scipy.special import digamma\n",
    "for volume in volumes:\n",
    "    for choice in range(count_samples):\n",
    "        x_sr = np.mean(data[volume][choice])\n",
    "        x_ln_sr = np.mean(np.log(data[volume][choice]))\n",
    "        eq = lambda theta: x_ln_sr + np.log(theta) - np.log(x_sr) - digamma(theta)\n",
    "        theta_1 = bisect(eq, 1, 10)\n",
    "        theta_2 = theta_1 / x_sr\n",
    "        print(f\"{'n = '+str(volume):>10},\",\n",
    "              f\"Выборка №{choice+1}, первый параметр = {theta_1:.3},\"\n",
    "              f\"второй параметр = {theta_2:.3}\")\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видна явная состоятельность оценки"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Продолжим исследование гамма распределения.<br>\n",
    "Теперь будем искать эффективную оценку для первого параметра распределения при известном втором.<br>\n",
    "$x_1, x_2, ... , x_n \\sim \\Gamma(\\theta, \\lambda)$, $\\theta$ - неизвестный параметр <br>\n",
    "$$V(\\vec{X}, \\theta)=\\frac{\\partial ln L}{\\partial \\theta} =  \\sum_{i=1}^{n}ln(x_i) + n \\space ln(\\lambda) - n \\space \\psi(\\theta) \\space \\space \\Leftrightarrow \\space\\space  \\frac{1}{n}V(\\vec{X}, \\theta) = \\frac{1}{n}\\sum_{i=1}^{n}ln(x_i) + ln(\\lambda) - \\psi(\\theta)$$ \n",
    "Таким образом эффективной оценкой для $\\psi(\\theta)$ является статистика $T(\\vec{X}) = \\frac{1}{n}\\sum_{i=1}^{n}ln(x_i) + ln(\\lambda) $, которая является фунцией от достаточной статистики, то есть оценка $\\widehat{\\psi(\\theta)}$ будет отпимальной.\n",
    "Обратная диагамма функция находится только численными методами."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Наконец будем искать эффективную оценку для второго параметра распределения при известном первом.<br>\n",
    "$x_1, x_2, ... , x_n \\sim \\Gamma(\\lambda, \\theta)$, $\\theta$ - неизвестный параметр <br>\n",
    "$$V(\\vec{X}, \\theta)=\\frac{\\partial ln L}{\\partial \\theta} =   n \\lambda \\frac{1}{\\theta} - \\sum_{i=1}^{n}x_i \\space \\space \\Leftrightarrow \\space\\space  \\frac{1}{n \\lambda}V(\\vec{X}, \\theta) = \\frac{1}{\\theta} - \\frac{1}{n \\lambda} \\sum_{i=1}^{n}x_i$$ \n",
    "Таким образом эффективной оценкой для $\\psi(\\theta)$ является статистика $T(\\vec{X}) = \\frac{1}{n \\lambda} \\sum_{i=1}^{n}x_i $, которая является фунцией от достаточной статистики, то есть оценка $\\frac{1}{\\hat{\\theta}}$ будет отпимальной.<br>\n",
    "Найдём оптимальную несмещённую оценку для $\\hat{\\theta}$. Для этого найдём смещение для $T_0(\\vec{X})= \\frac{n \\lambda}{\\sum_{i=1}^{n}x_i }$. Воспмним, что $\\sum_{i=1}^{n}x_i \\sim \\Gamma(n\\lambda, \\theta) $ <br>\n",
    "$$E[T_0] = \\int_{0}^{+\\infty} \\frac{n\\lambda}{x} \\frac{x^{n\\lambda-1}\\theta^{n\\lambda} e^{-x\\theta}}{\\Gamma(n\\lambda)} dx = \\frac{n\\lambda \\theta^{n\\lambda}}{\\Gamma(n\\lambda)} \\int_{0}^{+\\infty} (\\frac{z}{\\theta})^{n\\lambda-2} e^{-z} \\frac{dz}{\\theta} = \\theta  \\frac{n\\lambda}{\\Gamma(n\\lambda)} \\Gamma(n\\lambda - 1) = \\theta  \\frac{n\\lambda}{n\\lambda-1}$$\n",
    "Тогда несмещённой оптимальной оценкой будет $T(\\vec{X})=\\frac{n\\lambda-1}{n\\lambda}  \\frac{n \\lambda}{\\sum_{i=1}^{n}x_i } =  \\frac{n \\lambda -1 }{\\sum_{i=1}^{n}x_i } $ <br>\n",
    "Посмотрим, как ведёт себя эта оценка на смоделированных данных"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     n = 5, Выборка №1, искомый параметр = 0.498\n",
      "     n = 5, Выборка №2, искомый параметр = 0.631\n",
      "     n = 5, Выборка №3, искомый параметр = 0.407\n",
      "     n = 5, Выборка №4, искомый параметр = 0.393\n",
      "     n = 5, Выборка №5, искомый параметр = 0.691\n",
      "\n",
      "    n = 10, Выборка №1, искомый параметр = 1.21\n",
      "    n = 10, Выборка №2, искомый параметр = 0.449\n",
      "    n = 10, Выборка №3, искомый параметр = 0.409\n",
      "    n = 10, Выборка №4, искомый параметр = 0.572\n",
      "    n = 10, Выборка №5, искомый параметр = 0.727\n",
      "\n",
      "   n = 100, Выборка №1, искомый параметр = 0.6\n",
      "   n = 100, Выборка №2, искомый параметр = 0.435\n",
      "   n = 100, Выборка №3, искомый параметр = 0.449\n",
      "   n = 100, Выборка №4, искомый параметр = 0.548\n",
      "   n = 100, Выборка №5, искомый параметр = 0.528\n",
      "\n",
      "  n = 1000, Выборка №1, искомый параметр = 0.512\n",
      "  n = 1000, Выборка №2, искомый параметр = 0.499\n",
      "  n = 1000, Выборка №3, искомый параметр = 0.493\n",
      "  n = 1000, Выборка №4, искомый параметр = 0.49\n",
      "  n = 1000, Выборка №5, искомый параметр = 0.511\n",
      "\n",
      "n = 100000, Выборка №1, искомый параметр = 0.5\n",
      "n = 100000, Выборка №2, искомый параметр = 0.501\n",
      "n = 100000, Выборка №3, искомый параметр = 0.498\n",
      "n = 100000, Выборка №4, искомый параметр = 0.497\n",
      "n = 100000, Выборка №5, искомый параметр = 0.501\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for volume in volumes:\n",
    "    for choice in range(count_samples):\n",
    "        x_s = np.sum(data[volume][choice])\n",
    "        theta =  (volume*a - 1) / x_s\n",
    "        print(f\"{'n = '+str(volume):>10},\",\n",
    "              f\"Выборка №{choice+1}, искомый параметр = {theta:.3}\")\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видно, что даже при малых выборках значение близко к истинному"
   ]
  }
 ],
 "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.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}