ДЗ1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Гипергеометрическое распределение\n",
    "\n",
    "### Основные характеристики\n",
    "\n",
    "$\\xi \\sim HG (N, M, n), $ где $N$ - общее число элементов, $M$ - общее число элементов с нужным свойством, $n$ - размер выборки из общего числа элементов $ (n \\leq M \\leq N) $\n",
    "\n",
    "$P (\\xi = k) = \\frac{C_M^k C_{N-M}^{n-k}}{C_N^n}, k=0, ..., n $\n",
    "\n",
    "$ E[\\xi] = \\sum_{k=0}^{n}k\\frac{C_M^k C_{N-M}^{n-k}}{C_N^n} = \\sum_{k=1}^{n}k\\frac{M!C_{N-M}^{n-k}}{(M-k)!k!C_N^n} =  \\sum_{k=0}^{n-1}(k+1)\\frac{M!C_{N-M}^{n-k-1}}{(M-k-1)!(k+1)!C_N^n} = \\sum_{k=0}^{n-1}\\frac{M!C_{N-M}^{n-1-k}}{(M-1-k)!k!C_N^n} = M \\sum_{k=0}^{n-1}\\frac{C_{M-1}^k C_{N-M}^{n-1-k}}{C_N^n}  = $\n",
    "$= \\frac{M n}{N}\\sum_{k=0}^{n-1}\\frac{C_{M-1}^k C_{N-1-M+1}^{n-1-k}}{C_{N-1}^{n-1}} = \\frac{M n}{N}$\n",
    "\n",
    "Стоит пояснить, что $\\sum_{k=0}^{n-1}\\frac{C_{M-1}^k C_{N-1-M+1}^{n-1-k}}{C_{N-1}^{n-1}} = \\sum_{k=0}^{n-1}P(\\eta=k)=1,$ где $\\eta \\sim HG (N-1, M-1, n-1)$ \n",
    "\n",
    "\n",
    "$ E[\\xi^2] = \\sum_{k=0}^{n}k^2\\frac{C_M^k C_{N-M}^{n-k}}{C_N^n} = \\sum_{k=1}^{n}k^2\\frac{M!C_{N-M}^{n-k}}{(M-k)!k!C_N^n} =  \\sum_{k=0}^{n-1}(k+1)^2\\frac{M!C_{N-M}^{n-k-1}}{(M-k-1)!(k+1)!C_N^n} = \\sum_{k=0}^{n-1}(k+1)\\frac{M!C_{N-M}^{n-1-k}}{(M-1-k)!k!C_N^n} = $\n",
    "$ = M \\sum_{k=0}^{n-1}(k+1)\\frac{C_{M-1}^k C_{N-M}^{n-1-k}}{C_N^n}  = \\frac{M n}{N}\\sum_{k=0}^{n-1}(k+1)\\frac{C_{M-1}^k C_{N-1-M+1}^{n-1-k}}{C_{N-1}^{n-1}} = \\frac{M n}{N} E[\\eta + 1] = \\frac{M n}{N} (\\frac{(M-1)(n-1)}{(N-1)} + 1),$\n",
    "\n",
    "где $\\eta \\sim HG (N-1, M-1, n-1)$\n",
    "\n",
    "$D[\\xi] = E[\\xi^2] - (E[\\xi])^2 = \\frac{M n}{N} (\\frac{(M-1)(n-1)}{(N-1)} + 1 - \\frac{M n}{N})$\n",
    "\n",
    "Производящая функиця $g_\\xi(z) = E[z^\\xi] = \\sum_{k=0}^{n}z^k\\frac{C_M^k C_{N-M}^{n-k}}{C_N^n} $\n",
    "\n",
    "Характеристическая функция $\\phi_\\xi(t) = E[e^{it\\xi}] = \\sum_{k=0}^{n}e^{itk}\\frac{C_M^k C_{N-M}^{n-k}}{C_N^n} $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Зададим гипергеометрическое распределение $HG(200, 70, 10)$ используя модуль Scipy [1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import hypergeom\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "N, M, n = (200, 70, 10)\n",
    "hg = hypergeom(N, n, M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Построение эмпирических функций распределения"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAARYklEQVR4nO3dXYwd513H8e/PDmnZVkAge0Mce10IUIOggdO0UFFeWloXUMxFEQ5bFKCSBSK8SyjFF5WMLPEmKBcBakERgoVQAkIWooSKBu4avG5LwQmhrrGdxUVdcGmBRQ2u/1ycY3m9XWdns8dn9jz7/Uirc+aZmTP/0a5+OzPnmWdSVUiS2rWr7wIkSbeWQS9JjTPoJalxBr0kNc6gl6TG3dZ3AWvdeeedNTc313cZkjRVTp8+/e9VNbvevG0X9HNzcywuLvZdhiRNlSQXbjbPSzeS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6DU1FhZgbg527Rq+Liz0XZE0HbZd90ppPQsLcOQIrKwMpy9cGE4DzM/3V5c0DTyi11Q4evR6yF+zsjJsl/T8DHpNhYsXN9cu6TqDXlNh797NtUu6zqDXVDh+HGZmbmybmRm2S3p+nYI+ycEkzyQ5m+Thdeb/VJKnknw4yV8n2bdq3meSfGj0c3KcxWvnmJ+HEycAzgNX2bdvOO0XsdLGstEzY5PsBv4Z+DZgCTgFPFBVT61a5luAJ6tqJckPA99cVd8zmvffVfXSrgUNBoNyUDPdTBIAfNaxdKMkp6tqsN68Lkf09wFnq+pcVT0HPAocWr1AVT1RVdf6RLwf2LOVgiVJ49Ml6O8Cnl01vTRqu5m3Au9ZNf3iJItJ3p/ku9ZbIcmR0TKLy8vLHUqSJHXV5YaprNO27nlzkrcAA+CbVjXvrapLSV4GvC/JP1TVR2/4sKoTwAkYXrrpVLkkqZMuR/RLwN2rpvcAl9YulOT1wFHg/qr69LX2qro0ej0H/A1w7xbqlXrh8AuaZl2C/hRwT5L9SW4HDgM39J5Jci/wToYh//FV7XckedHo/Z3Aa4CnkKbIteEXLlyAquvDLxj2mhYbBn1VXQEeAh4HngbeXVVnkhxLcv9osV8CXgr88ZpulC8HFpP8PfAE8POre+tI08DhFzTtNuxeOWl2r9Tz6aN75a5dwyP5z64Frl6dWBnS89pq90ppR3P4BU07g17agMMvaNoZ9NIGHH5B085r9JoqfQ6B4PAL2s68Ri9JO5hBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJalynoE9yMMkzSc4meXid+T+V5KkkH07y10n2rZr3YJKPjH4eHGfxkqSNbRj0SXYDjwBvAg4ADyQ5sGaxDwKDqvpq4DHgF0frfiHwduBVwH3A25PcMb7yJUkb6XJEfx9wtqrOVdVzwKPAodULVNUTVbUymnw/sGf0/o3Ae6vqclV9AngvcHA8pUuSuugS9HcBz66aXhq13cxbgfdsZt0kR5IsJllcXl7uUJIkqasuQZ912mrdBZO3AAPglzazblWdqKpBVQ1mZ2c7lCRJ6qpL0C8Bd6+a3gNcWrtQktcDR4H7q+rTm1lXknTrdAn6U8A9SfYnuR04DJxcvUCSe4F3Mgz5j6+a9TjwhiR3jL6EfcOoTZI0IbdttEBVXUnyEMOA3g28q6rOJDkGLFbVSYaXal4K/HESgItVdX9VXU7ycwz/WQAcq6rLt2RPJEnrStW6l9t7MxgManFxse8ytE2NDiTo4++2z21LG0lyuqoG683zzlhJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr02bWEB5uZg167h68JC3xVJej4bjl4prbawAEeOwMrowZEXLgynAebn+6tL0s15RK9NOXr0eshfs7IybJe0PRn02pSLFzfXLql/Br02Ze/ezbVL6p9Br005fhxmZm5sm5kZtkvangx6bcr8PJw4AXAeuMq+fcNpv4iVti8fJagXpK/H6vkoQWl9PkpQknYwg16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXtjGHhNY4OEyxtE05JLTGxSN6aZtySGiNi0EvbVMOCa1x6RT0SQ4meSbJ2SQPrzP/tUk+kORKkjevmfeZJB8a/ZwcV+FS6xwSWuOyYdAn2Q08ArwJOAA8kOTAmsUuAt8P/ME6H/G/VfWK0c/9W6xX2jEcElrj0uWI/j7gbFWdq6rngEeBQ6sXqKrzVfVh4OotqFHakRwSWuPSJejvAp5dNb00auvqxUkWk7w/yXett0CSI6NlFpeXlzfx0VLbhqG+H9jN+fOGvF6YLkGfddo2MyD33tEYyd8LvCPJl3zWh1WdqKpBVQ1mZ2c38dGSpI10Cfol4O5V03uAS103UFWXRq/ngL8B7t1EfZKkLeoS9KeAe5LsT3I7cBjo1HsmyR1JXjR6fyfwGuCpF1qsJGnzNgz6qroCPAQ8DjwNvLuqziQ5luR+gCSvTLIEfDfwziRnRqu/HFhM8vfAE8DPV5VBL0kT5DNj9YL4zNj2t6vp4jNjJWkHM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUuE5Bn+RgkmeSnE3y8DrzX5vkA0muJHnzmnkPJvnI6OfBcRUuSepmw6BPsht4BHgTcAB4IMmBNYtdBL4f+IM1634h8HbgVcB9wNuT3LH1siVJXXU5or8POFtV56rqOeBR4NDqBarqfFV9GLi6Zt03Au+tqstV9QngvcDBMdQtSeqoS9DfBTy7anpp1NZFp3WTHEmymGRxeXm540dLkrroEvRZp606fn6ndavqRFUNqmowOzvb8aMlSV10Cfol4O5V03uASx0/fyvrSpLGoEvQnwLuSbI/ye3AYeBkx89/HHhDkjtGX8K+YdQmSZqQDYO+qq4ADzEM6KeBd1fVmSTHktwPkOSVSZaA7wbemeTMaN3LwM8x/GdxCjg2apMkTUiqul5un4zBYFCLi4t9l6ENJMOvXyb999PXdvvcdp/7rOmR5HRVDdab552xktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gn2ILCzA3B7t2DV8XFvquSNJ2dFvfBeiFWViAI0dgZWU4feHCcBpgfr6/uiRtPx7RT6mjR6+H/DUrK8N2aRw8Y2yHR/RT6uLFzbVLm+EZY1s8op9Se/durl3aDM8Y22LQT6njx2Fm5sa2mZlhu7RVnjG2xaCfUvPzcOIEwHngKvv2Dac9rdY4eMbYFoN+ig1DfT+wm/PnDXmNj2eMbTHoJX0Wzxjbkqrqu4YbDAaDWlxc7LuMqZEEgEn/Hnfadvvc9k7cZ21ektNVNVhvnkf0ktQ4g16SGmfQS1LjDHpJalynoE9yMMkzSc4meXid+S9K8kej+U8mmRu1zyX53yQfGv385njLlyRtZMOxbpLsBh4Bvg1YAk4lOVlVT61a7K3AJ6rqS5McBn4B+J7RvI9W1SvGXLckqaMuR/T3AWer6lxVPQc8Chxas8wh4HdH7x8DXpdr/bIkSb3qEvR3Ac+uml4ata27TFVdAT4JfNFo3v4kH0zyt0m+cYv1SpI2qcswxesdma+9e+Jmy3wM2FtV/5Hk64A/S/KVVfWpG1ZOjgBHAPY6mIYkjVWXI/ol4O5V03uASzdbJsltwOcDl6vq01X1HwBVdRr4KPBlazdQVSeqalBVg9nZ2c3vhSTpproE/SngniT7k9wOHAZOrlnmJPDg6P2bgfdVVSWZHX2ZS5KXAfcA58ZTuiSpiw0v3VTVlSQPAY8Du4F3VdWZJMeAxao6Cfw28HtJzgKXGf4zAHgtcCzJFeAzwA9V1eVbsSOSpPU5qNmU22kDbe3EAb524j5r8xzUTJJ2MINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJe0rSwswNwc7No1fF1Y6Lui6dfl4eCSNBELC3DkCKysDKcvXBhOA8zP91fXtPOIXtK2cfTo9ZC/ZmVl2K4XzqCXtG1cvLi5dnVj0EvaNvbu3Vy7ujHoJW0bx4/DzMyNbTMzw3a9cAb9FtlDQBqf+Xk4cQLgPHCVffuG034RuzWpqr5ruMFgMKjFxcW+y+hkbQ8BGB59TPIPMwkAk/497rTt9rlt91ldJDldVYP15nlEvwX2EJA0DQz6LbCHgKRpYNBvgT0EJE0Dg34L7CEgaRoY9FtgDwGpLa32orPXzRjYK6L97fa5bfd5MrZDL7qtsNeNJG2g5V50nYI+ycEkzyQ5m+Thdea/KMkfjeY/mWRu1by3jdqfSfLG8ZV+o1ZPuSRNRp+96G51fm0Y9El2A48AbwIOAA8kObBmsbcCn6iqLwV+FfiF0boHgMPAVwIHgV8ffd5YXTvlunABqq4PbWrYS+qqr150k8ivLkf09wFnq+pcVT0HPAocWrPMIeB3R+8fA16X4UW2Q8CjVfXpqvoX4Ozo88aq5VMuSZPRVy+6SeRXl6C/C3h21fTSqG3dZarqCvBJ4Is6rkuSI0kWkywuLy93r36k7xuXqqq3W7X72vZO226f23afJ+NaL7p9+yBhYr3oJpFfXYI+67St/Q3cbJku61JVJ6pqUFWD2dnZDiXdyBuXJI3D/DycPw9Xrw5fJ9HbZhL51SXol4C7V03vAS7dbJkktwGfD1zuuO6WeeOSpGk1ifzqEvSngHuS7E9yO8MvV0+uWeYk8ODo/ZuB99XwvOskcHjUK2c/cA/wd+Mp/bq+TrkkaasmkV+dbphK8u3AO4DdwLuq6niSY8BiVZ1M8mLg94B7GR7JH66qc6N1jwI/CFwBfqKq3vN825rGG6YkqW/Pd8OUd8ZKUgO8M1aSdjCDXpIaZ9BLUuMMeklq3Lb7MjbJMnCh7zpegDuBf++7iAlzn3cG93k67Kuqde843XZBP62SLN7sG+9Wuc87g/s8/bx0I0mNM+glqXEG/fic6LuAHrjPO4P7POW8Ri9JjfOIXpIaZ9BLUuMM+jHY6OHprUlyd5Inkjyd5EySH++7pklIsjvJB5P8ed+1TEKSL0jyWJJ/Gv2uv77vmm61JD85+pv+xyR/OBqZd+oZ9FvU8eHprbkC/HRVvRx4NfAjO2CfAX4ceLrvIibo14C/rKqvAL6Gxvc9yV3AjwGDqvoqhsOyH+63qvEw6Leuy8PTm1JVH6uqD4ze/xfDAPisZwG3JMke4DuA3+q7lklI8nnAa4HfBqiq56rqP/utaiJuAz539KS8GW7BE/H6YNBvXacHoLcqyRzDB8482W8lt9w7gJ8BrvZdyIS8DFgGfmd0ueq3kryk76Jupar6V+CXgYvAx4BPVtVf9VvVeBj0W9fpAegtSvJS4E8YPjnsU33Xc6sk+U7g41V1uu9aJug24GuB36iqe4H/AZr+/inJHQzPxvcDXwy8JMlb+q1qPAz6rZvIA9C3mySfwzDkF6rqT/uu5xZ7DXB/kvMML819a5Lf77ekW24JWKqqa2dqjzEM/pa9HviXqlquqv8D/hT4hp5rGguDfuu6PDy9KUnC8Nrt01X1K33Xc6tV1duqak9VzTH8/b6vqpo40ruZqvo34NkkXz5qeh3wVI8lTcJF4NVJZkZ/46+jkS+gb+u7gGlXVVeSPAQ8zvWHp5/puaxb7TXA9wH/kORDo7afraq/6LEmjd+PAgujA5hzwA/0XM8tVVVPJnkM+ADDnmUfpJGhEBwCQZIa56UbSWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIa9//uVAyxhhQpnwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.arange(-1, n)\n",
    "pmf = hg.pmf(x)\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(x, pmf, 'bo')\n",
    "ax.vlines(x, 0, pmf, lw=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Функция распределения"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAR80lEQVR4nO3df6yddX3A8feH1pbxQ9nsZWpvSzGrw2pE3IWxkSgIxqKjXYzTNnE6Q2gIojjJFtQFCYuJgznnMtRVhqhTCGNmNqZbtyiLZhnY8kNm26FNkfYONq4Oma5gpf3sj3OuO7u9t/dw+zzPOef7vF9Jc8+Ph/N8T3rum2+/5znnicxEkjT6jhv0ACRJ1TDoklQIgy5JhTDoklQIgy5JhVg8qB0vW7YsV61aNajdS9JIuvfee7+fmWOz3TewoK9atYodO3YMaveSNJIi4pG57nPJRZIKYdAlqRAGXZIKYdAlqRAGXZIKMW/QI+KWiHg8Ir49x/0REX8WEXsi4sGIeFX1w5QkzaefGfqtwNqj3H8xsLr7ZxPwyWMfliTp2Zr3OPTM/HpErDrKJuuBz2Xne3jvjohTIuKFmflYRWOU1CI/+Qnceis89VTneu83fFd9uc7HPtrlSy6Bs8+mclV8sGg5sL/n+mT3tiOCHhGb6MziWblyZQW7ltSEq6+Gb32rmX3dcw/8+MfN7KtpEZ2fL3rR8AY9Zrlt1rNmZOZmYDPAxMSEZ9aQRkAmfOxjMD4OTczDzjwTxsbg05+Gxd1CRU9lqr5c52M3rYqgTwIreq6PA49W8LiShsDBg52oX345fOADgx6NjqaKwxa3AG/vHu1yLvCk6+dSOQ4c6Pw84YTBjkPzm3eGHhG3AecDyyJiEvgQ8ByAzPwUsBV4A7AHOAC8s67BSuq48Ua47rpm9jX9Zt6JJzazPy1cP0e5bJzn/gTeVdmIJM3rvvtg6VK49NJm9rdkCaxf38y+tHAD+/pcSQt36BC84AWdmbo0zY/+SyPo8GE4zt9ezeBLQhpBhw7BokWDHoWGjUGXRpBB12wMujSCDh1yyUVH8iUhjaDDh52h60gGXRpBLrloNgZdGkEuuWg2viSkEeQMXbMx6NIIcg1dszHo0ghyyUWz8aP/UgW++114y1v+7yw7dXvkEXj1q5vZl0aHQZcq8OCD8MAD8PrXwymn1L+/V74SNh71a/PURgZdqsD0V8zecAO84hWDHYvay1U4qQLTQR/UqcckMOhSJQy6hoFBlyowHXSPPNEg+fKTKuAMXcPAoEsVMOgaBgZdqsDhw52fBl2DZNClCjhD1zAw6FIFDLqGgUGXKuBRLhoGvvykCjhD1zAw6FIFDLqGgUGXKmDQNQwMulQBD1vUMDDoUgV8U1TDwJefVAGXXDQMDLpUAYOuYWDQpQoYdA2DvoIeEWsj4qGI2BMR18xy/8qIuCsi7o+IByPiDdUPVRpeBl3DYN6gR8Qi4CbgYmANsDEi1szY7A+AOzLzLGAD8ImqByoNM4OuYdDPDP0cYE9m7s3Mg8DtwPoZ2yTw3O7l5wGPVjdEafhNH7boUS4apH5efsuB/T3XJ7u39boOeFtETAJbgXfP9kARsSkidkTEjqmpqQUMVxpOztA1DPoJ+mwv0ZxxfSNwa2aOA28APh8RRzx2Zm7OzInMnBgbG3v2o5WGlEHXMOgn6JPAip7r4xy5pHIpcAdAZv4LcDywrIoBSqPAoGsY9BP07cDqiDg9IpbQedNzy4xt9gEXAkTES+kE3TUVtYZB1zCYN+iZ+QxwJbAN2E3naJadEXF9RKzrbnY1cFlEfAu4DfidzJy5LCMVy6BrGCzuZ6PM3Ernzc7e267tubwLOK/aoUmjw+9y0TDw5SdVwG9b1DDoa4YujaJvfhP27m1mX/ff3/lp0DVIBl3FuuACOHCguf2ddBIsXdrc/qSZDLqKdeAAXH45XHVVM/tbtgyWLGlmX9JsDLqKNjYGZ5wx6FFIzfBNURXNNW21iUFXkTwuXG1k0FUkP9amNjLoKpozdLWJQVeRXHJRGxl0Fcmgq40MuiQVwqCrSM7Q1UYGXUUy6Gojg64iGXS1kUGXpEIYdBXJGbrayKCrSAZdbWTQJakQBl1FcoauNjLoKpJBVxsZdBXJoKuNDLokFcKgq0jO0NVGBl1FMuhqI4Ouohl0tYlBV5E8BZ3ayKCrSC65qI0Muopk0NVGBl1FM+hqk76CHhFrI+KhiNgTEdfMsc1bImJXROyMiC9WO0zp2XENXW20eL4NImIRcBPwOmAS2B4RWzJzV882q4H3A+dl5hMRcWpdA5b64ZKL2qifGfo5wJ7M3JuZB4HbgfUztrkMuCkznwDIzMerHaa0MAZdbdJP0JcD+3uuT3Zv6/US4CUR8c8RcXdErJ3tgSJiU0TsiIgdU1NTCxux1Adn6GqjfoI+26/EzBXKxcBq4HxgI3BzRJxyxH+UuTkzJzJzYmxs7NmOVeqba+hqo36CPgms6Lk+Djw6yzZfzsyfZubDwEN0Ai8NhDN0tVE/Qd8OrI6I0yNiCbAB2DJjm78FLgCIiGV0lmD2VjlQaSEMutpk3qBn5jPAlcA2YDdwR2bujIjrI2Jdd7NtwA8iYhdwF/B7mfmDugYtzcclF7XRvIctAmTmVmDrjNuu7bmcwPu6f6SBc8lFbeQnRVU0g642MegqkjN0tZFBV5FcQ1cbGXQVyRm62sigq2gGXW1i0FUkZ+hqI4OuIrmGrjYy6CqaM3S1iUFXkVxyURsZdBXJoKuN+vrov1SFp56CXbvm364Kjz3WzH6kYWLQ1Zj3vhc2b252nyee2Oz+pEEy6GrME0/A+Dh84hPN7G/pUrjggmb2JQ0Dg67GZMLJJ8Mllwx6JFKZfFNUjcn0TUqpTgZdjfHDPlK9DLoa5Qxdqo9BV2NccpHqZdDVGJdcpHoZdDXKGbpUH4OuxrjkItXLoKsxBl2ql0GXpEIYdDXGGbpUL4Ouxhh0qV4GXZIKYdDVGGfoUr0Muhpj0KV6GXRJKoRBV2OcoUv1MuhqjEGX6tVX0CNibUQ8FBF7IuKao2z35ojIiJiobogqiUGX6jNv0CNiEXATcDGwBtgYEWtm2e5k4D3APVUPUmXw2xalevUzQz8H2JOZezPzIHA7sH6W7f4QuAF4usLxqSAuuUj16ifoy4H9Pdcnu7f9TEScBazIzK8c7YEiYlNE7IiIHVNTU896sBp9Bl2qTz9Bn+1X8Gf/eI6I44CPAVfP90CZuTkzJzJzYmxsrP9RqgguuUj16ifok8CKnuvjwKM9108GXg78U0R8DzgX2OIbo5rJJRepXv0EfTuwOiJOj4glwAZgy/SdmflkZi7LzFWZuQq4G1iXmTtqGbFGmkGX6jNv0DPzGeBKYBuwG7gjM3dGxPURsa7uAaocLrlI9Vrcz0aZuRXYOuO2a+fY9vxjH5ZK5JKLVC8/KapGGXSpPgZdjXGGLtXLoKsxrqFL9TLoapQzdKk+Bl2NcclFqpdBV2NccpHqZdDVGGfoUr0Muhpl0KX6GHQ1xiUXqV4GXY1xyUWql0FXowy6VB+Drsa45CLVy6CrMS65SPUy6GqUQZfqY9DVGGfoUr0MuhrjGrpUL4OuRjlDl+pj0NUYl1ykehl0NcYlF6leBl2NcoYu1cegqzEuuUj1WjzoAWhwDh2Cb3wDnnqqmf398Idw6qnN7EtqI4PeYtu2wRvf2Ow+zz672f1JbWLQW+xHP+r8/OIX4cUvbmafL3tZM/uR2sigt9j0USdnnglr1gx2LJKOnW+Ktth00H2jUiqDQW8xgy6VxaC3mEGXymLQW8ygS2Ux6C02HfTjfBVIRejrVzki1kbEQxGxJyKumeX+90XEroh4MCK+GhGnVT9UVc0ZulSWeYMeEYuAm4CLgTXAxoiYeZDb/cBEZr4CuBO4oeqBqnoGXSpLPzP0c4A9mbk3Mw8CtwPrezfIzLsy80D36t3AeLXDVB0MulSWfoK+HNjfc32ye9tcLgX+brY7ImJTROyIiB1TU1P9j1K1MOhSWfoJ+my/7rN+s3VEvA2YAG6c7f7M3JyZE5k5MTY21v8oVQuDLpWln4/+TwIreq6PA4/O3CgiLgI+CLwmM39SzfBUJ4MulaWfGfp2YHVEnB4RS4ANwJbeDSLiLOAvgHWZ+Xj1w1QdDLpUlnmDnpnPAFcC24DdwB2ZuTMiro+Idd3NbgROAv46Ih6IiC1zPJyGiEGXytLXty1m5lZg64zbru25fFHF41IDDLpUFj8j2GIGXSqLQW8xgy6VxaC3mEGXymLQW8ygS2Ux6C1m0KWyGPQWM+hSWQx6ixl0qSwGvcUMulQWg95iBl0qi0FvMYMulcWgt5hBl8pi0FvMoEtlMegtZtClshj0FjPoUlkMeosZdKksBr3FDLpUFoPeYgZdKktfZyxSMzLhiitg9+5m9rdvX+enQZfKYNCHyOHD8KlPwWmnwapV9e9v5Up47Wvh+OPr35ek+hn0IXLoUOfnZZfBBz842LFIGj2uoQ+Rw4c7PxctGuw4JI0mgz5Epmfox/m3ImkBTMcQcYYu6VgY9CHiDF3SsTAdQ2Q66M7QJS2EQR8iLrlIOhYGfYi45CLpWJiOIeIMXdKxMOhDxBm6pGNhOoaIM3RJx8KP/s/j6afhO99pZl+Tk52fztAlLURfQY+ItcDHgUXAzZn5kRn3LwU+B/wK8APgrZn5vWqHOhhXXAGf+Uyz+zzppGb3J6kM8wY9IhYBNwGvAyaB7RGxJTN39Wx2KfBEZv5SRGwA/gh4ax0Dbtq+fXDGGfDhDzezv+OPh4suamZfksrSzwz9HGBPZu4FiIjbgfVAb9DXA9d1L98J/HlEROb0KRSqc8st8NGPVv2oc3v4YXjNa+BNb2pun5K0EP0EfTmwv+f6JPCrc22Tmc9ExJPA84Hv924UEZuATQArV65c0ICf/3xYs2ZB/+mCrFkDb397c/uTpIXqJ+iznc9m5sy7n23IzM3AZoCJiYkFzd7Xr+/8kST9f/0cTzEJrOi5Pg48Otc2EbEYeB7wX1UMUJLUn36Cvh1YHRGnR8QSYAOwZcY2W4B3dC+/GfhaHevnkqS5zbvk0l0TvxLYRuewxVsyc2dEXA/syMwtwF8Cn4+IPXRm5hvqHLQk6Uh9HYeemVuBrTNuu7bn8tPAb1U7NEnSs+FnEiWpEAZdkgph0CWpEAZdkgoRgzq6MCKmgEcGsvNnbxkzPvVakJKfG5T9/Hxuo+tYnt9pmTk22x0DC/ooiYgdmTkx6HHUoeTnBmU/P5/b6Krr+bnkIkmFMOiSVAiD3p/Ngx5AjUp+blD28/O5ja5anp9r6JJUCGfoklQIgy5JhTDoRxERayPioYjYExHXDHo8VYqIFRFxV0TsjoidEXHVoMdUtYhYFBH3R8RXBj2WqkXEKRFxZ0T8W/fv8NcGPaaqRMTvdl+T346I2yLi+EGPaaEi4paIeDwivt1z2y9ExD9GxHe7P3++qv0Z9Dn0nBz7YmANsDEiGjz5Xe2eAa7OzJcC5wLvKuz5AVwF7B70IGryceDvM/MM4EwKeZ4RsRx4DzCRmS+n85Xdo/x13LcCa2fcdg3w1cxcDXy1e70SBn1uPzs5dmYeBKZPjl2EzHwsM+/rXv4RnSAsH+yoqhMR48AbgZsHPZaqRcRzgVfTOQ8BmXkwM3842FFVajHwc92zn53AkWdIGxmZ+XWOPHvbeuCz3cufBX6zqv0Z9LnNdnLsYoLXKyJWAWcB9wx2JJX6U+D3gcODHkgNXgxMAZ/pLindHBEnDnpQVcjMfwf+GNgHPAY8mZn/MNhRVe4XM/Mx6EysgFOremCDPre+Tnw96iLiJOBvgPdm5n8PejxViIjfAB7PzHsHPZaaLAZeBXwyM88C/ocK/9k+SN315PXA6cCLgBMj4m2DHdXoMOhz6+fk2CMtIp5DJ+ZfyMwvDXo8FToPWBcR36OzVPbaiPirwQ6pUpPAZGZO/4vqTjqBL8FFwMOZOZWZPwW+BPz6gMdUtf+MiBcCdH8+XtUDG/S59XNy7JEVEUFnDXZ3Zv7JoMdTpcx8f2aOZ+YqOn9vX8vMYmZ5mfkfwP6I+OXuTRcCuwY4pCrtA86NiBO6r9ELKeQN3x5bgHd0L78D+HJVD9zXOUXbaK6TYw94WFU6D/ht4F8j4oHubR/onj9Ww+/dwBe6k429wDsHPJ5KZOY9EXEncB+dI7HuZ4S/BiAibgPOB5ZFxCTwIeAjwB0RcSmd/4FVdj5mP/ovSYVwyUWSCmHQJakQBl2SCmHQJakQBl2SCmHQJakQBl2SCvG/1Fbtr+FGfvgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.arange(-1, n, 0.001)\n",
    "cdf = hg.cdf(x)\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(x, cdf, 'b')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Классической моделью применения является оценка численности замкнутой популяции животных.\n",
    "\n",
    "    Проблема состоит, например, в определении численности N рыб, живущих на момент наблюдения в замкнутом водоеме, скажем, в пруду рыбоводного хозяйства [2]. Для определения (точнее, приближенной оценки) N исследователь отлавливает заданное количество M рыб, метит их каким-либо способом и возвращает в пруд. По истечении некоторого промежутка времени, когда, по его мнению, меченые рыбы “перемешались” с другими обитателями пруда, он снова отлавливает фиксированное количество n рыб (в математической статистике эта процедура называется извлечением выборки объема n из генеральной совокупности) и подсчитывает число m отмеченных рыб, попавших во второй улов. В рамках гипергеометрической модели такого эксперимента мы располагаем значениями параметров M и n, знаем результат m наблюдения случайной величины X, но не знаем значения параметра N гипергеометрического распределения P(X = k | N, M, n).\n",
    "\n",
    "    Один из основных методов решения обратных задач теории вероятностей (задач математической статистики), который называется методом максимального правдоподобия, состоит в выборе такого значения Nˆ параметра N, которое соответствует максимуму вероятности наблюдаемого исхода m в наблюдении X. Основной довод в пользу такого поведения статистика состоит в простом житейском наблюдении: если происходит какоелибо событие, то это событие должно иметь большую вероятность по сравнению с другими исходами статистического эксперимента.\n",
    "\n",
    "    Итак, метод максимального правдоподобия предлагает в качестве оценки неизвестного значения N (численности рыб в пруду) взять решения следующей задачи на экстремум:\n",
    "\n",
    "<center>$\\hat{N} = arg \\underset{N}\\max{\\frac{C_M^k C_{N-M}^{n-k}}{C_N^n}}  $</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    Решить эту задачу можно с помощью определения значения N, при котором происходит смена неравенства\n",
    "<center>$\\hat{N} = \\frac{C_M^k C_{N-M}^{n-k}}{C_N^n} < \\frac{C_M^k C_{N+1-M}^{n-k}}{C_{N+1}^n}  $</center>\n",
    "    \n",
    "    на обратное. Используя известную формулу для вычисления биномиальных коэффициентов, находим, что это неравенство эквивалентно (N + 1)k > nM, откуда получаем оценку максимального правдоподобия для численности рыб в пруду:\n",
    "\n",
    "<center>$\\hat{N} = k   \\lfloor \\frac{n}{M} \\rfloor $</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Моделирование гипергеометрического распределения"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Вариантом моделирования случайно величины $\\xi$ с известной функцией распределения $F_\\xi(x)$ является метод обратных функций.[3]\n",
    "\n",
    "Заметим, что $F_\\xi: \\mathbb{R} \\rightarrow [0,1]$\n",
    "\n",
    "Пусть $\\eta \\sim R[0,1]$ - случайная величина, которую мы умеем моделировать, а $F_\\eta (x) = x$ - её функция распределения.\n",
    "\n",
    "Тогда $ F_{F_\\xi^{-1}(\\eta)}(x) = P(F_\\xi^{-1}(\\eta)<x) = P(\\eta<F_\\xi(x)) = F_\\eta(F_\\xi(x)) = F_\\xi(x)$\n",
    "\n",
    "Иными словами, если у нас есть $X_1, ..., X_n$ - выборка из стандартного равномерного распределения $R[0,1]$, то \n",
    "\n",
    "$F_\\xi^{-1} (X_1), ..., F_\\xi^{-1} (X_n) $ - выборка из моделируемого распределения случайно величины $\\xi$.\n",
    "\n",
    "Стоит заметить, что утверждение выше наклдывает ограничение на функцию $F_\\xi(x)$, которая должна быть биективной (иными словами строговозрастающей). Функция распределения гипергеометрического распределения является неубывающей, но никак не строговозрастающей. Для того, чтобы решить этот вопрос будет брать точную нижнюю границу полного праоброза известной нам функции распределения, то есть $inf\\{x| F_\\xi(x) = X_1\\}, ... ,inf\\{x| F_\\xi(x) = X_n\\}$- выборка моделируемого распределения случайно величины $\\xi$ свободной от ограничений на функцию распределения. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Гамма распределение\n",
    "\n",
    "### Основные характеристики\n",
    "\n",
    "$\\xi \\sim \\Gamma(\\alpha, \\theta), $\n",
    "\n",
    "$f(x) = \\frac{x^{\\alpha-1}\\theta^\\alpha e^{-x\\theta}}{\\Gamma(\\alpha)}, x \\geq 0$\n",
    "\n",
    "$ E[\\xi] = \\int_{-\\infty }^{+\\infty}x f(x)dx = \\int_{0}^{+\\infty}x \\frac{x^{\\alpha-1}\\theta^\\alpha e^{-x\\theta}}{\\Gamma(\\alpha)} dx= \\frac{\\theta^\\alpha}{\\Gamma(\\alpha)} \\int_{0}^{+\\infty} x^{\\alpha}e^{-x\\theta} dx = \\frac{\\theta^\\alpha}{\\Gamma(\\alpha)} \\int_{0}^{+\\infty}  \\big(\\frac{z}{\\theta} \\big)^\\alpha e^{-z} \\frac{dz}{\\theta} = \\frac{1}{\\theta\\ \\Gamma(\\alpha)}\\int_{0}^{+\\infty}z^\\alpha e^{-z}dz=\\frac{\\Gamma(\\alpha+1)}{\\theta\\ \\Gamma(\\alpha)}=\\frac{\\alpha}{\\theta}$\n",
    "\n",
    "$ E[\\xi^2] = \\int_{-\\infty }^{+\\infty}x^2 f(x)dx = \\int_{0}^{+\\infty}x^2 \\frac{x^{\\alpha-1}\\theta^\\alpha e^{-x\\theta}}{\\Gamma(\\alpha)} dx= \\frac{\\theta^\\alpha}{\\Gamma(\\alpha)} \\int_{0}^{+\\infty} x^{\\alpha+1}e^{-x\\theta} dx = \\frac{\\theta^\\alpha}{\\Gamma(\\alpha)} \\int_{0}^{+\\infty}  \\big(\\frac{z}{\\theta} \\big)^\\alpha+1 e^{-z} \\frac{dz}{\\theta} = \\frac{1}{\\theta^2\\ \\Gamma(\\alpha)}\\int_{0}^{+\\infty}z^\\alpha+1 e^{-z}dz=\\frac{\\Gamma(\\alpha+2)}{\\theta^2\\ \\Gamma(\\alpha)}=\\frac{\\alpha(\\alpha+1)}{\\theta^2}$\n",
    "\n",
    "$D[\\xi] = E[\\xi^2] - (E[\\xi])^2 = \\frac{\\alpha}{\\theta^2}$\n",
    "\n",
    "Характеристическая функция $\\phi_\\xi(t) = E[e^{it\\xi}] =  \\int_{0}^{+\\infty} e^{itx} \\frac{x^{\\alpha-1}\\theta^\\alpha e^{-x\\theta}}{\\Gamma(\\alpha)} dx= \\frac{\\theta^\\alpha}{\\Gamma(\\alpha)} \\int_{0}^{+\\infty} e^{(it-\\theta) x} x^{\\alpha-1} dx= \\frac{\\theta^\\alpha}{\\Gamma(\\alpha)} \\int_{0}^{+\\infty} e^{z} (\\frac{z}{it-\\theta})^{\\alpha-1} \\frac{dz}{it-\\theta}= \\frac{\\theta^\\alpha (it-\\theta)^{-\\alpha} \\Gamma(\\alpha)}{\\Gamma(\\alpha)} = \\theta^\\alpha (it-\\theta)^{-\\alpha} $\n",
    "\n",
    "### Известные соотношения:\n",
    "\n",
    "Если $\\xi_1, ... , \\xi_n $ - независимые экспоненциальные случайные величины, такие что  $ \\xi_i \\sim Exp(\\theta) \\space\\forall i=1,..,n$,то $ \\sum_{i=1}^n \\xi_i \\sim \\Gamma(n, \\theta)$\n",
    "\n",
    "Если $ \\xi \\sim \\Gamma(\\alpha, \\theta)$, то $c \\xi  \\sim \\Gamma(\\alpha, c\\theta)  \\space\\forall c>0$\n",
    "\n",
    "$\\Gamma(n/2, 1/2) \\equiv \\chi^2 (n), n \\in \\mathbb{N} $\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Применения гамма-распределения"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Классическим методом применения гамма-распределения является моделирование других случайных величин."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Помимо этого гамма-распледелением описывает уязвимость живого оргазима к кумлятивному риску смерти. [4]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Гамма-распределению могут быть подчинены такие величины, как общий срок службы изделия, время наработки до k-го отказа.\n",
    "Также, это распределение используется в логистике для описания спроса в моделях управления запасами [5]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Выведем плотность гамма распределения $\\Gamma(2, 1/2)$, задав эти функции вручную"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy.integrate as integrate\n",
    "import scipy.special as sps\n",
    "%matplotlib inline\n",
    "a, th = (2, 1/2)\n",
    "def f_gamma(x):\n",
    "    return np.exp(-x*th) * th**a * x**(a-1) / (sps.gamma(a))\n",
    "\n",
    "def F_gamma(x):\n",
    "    return [(integrate.quad(f_gamma, 0, a)[0] if a>0 else 0) for a in x]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Плотность распределения"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD4CAYAAADlwTGnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de3xU9Zn48c8zM7knkIQkCCRAuBREEcSAgECtthYvFd1qxVpBq9Xa+mv319/u1nZb27Xb3W3dtttWbWu9WxVdrZVaLLVeWq9AUO4XDeEWCCGQ+/32/P6YExzHCZnATE6S87xfr3nNme/5nnOemUzOM9/vOed7RFUxxhjjPT63AzDGGOMOSwDGGONRlgCMMcajLAEYY4xHWQIwxhiPCrgdQF/k5OTo+PHj3Q7DGGMGlfXr1x9R1dzw8kGVAMaPH09xcbHbYRhjzKAiInsjlVsXkDHGeJQlAGOM8ShLAMYY41GWAIwxxqMsARhjjEdZAjDGGI+yBGCMMR41qK4DGEwO1jTzt/cqOVLfSm5GEudMyqEgO9XtsIwx5hhLADFW29zOf/xpO0+t30/4rRYWn3YK377oVMaOsERgjHGfJYAYOlTbwrIH1lBa2cgXzynk82ePJT8rhbLqZp579wAPvLGHi37xGj++4gwumj7K7XCNMR5nxwBipKG1g+seXMvBmhYe+eIcvnvJNCbmppMU8DMxN51vXDCFP//jQj42Mp2vPv4Oj74d8cpsY4zpN5YAYuQ7z27m/cMN3HPNLOZPyolYJz8rlce/NJfzp47ku3/YwtPry/o5SmOM+YAlgBh4aXsFf9hwkP9z3iQWfewjA+59SHKCn3uumcU5k0Zw2zObeP39I/0UpTHGfJglgJPU2tHJ7c9tZcrIDL5y7qSolkkM+PjVF85iYm46X1vxLuW1zXGO0hhjPsoSwEl6ct1+DtQ08+2LTyUxEP3HOSw5gXu+MIuW9k6+9sS7dHR2xTFKY4z5qKj2WCKyWER2ikiJiNwWYf4iEXlHRDpE5IqQ8k+IyIaQR4uIXObMe0hEdofMmxm7t9U/Wjs6uevlEuaMz2bR5Mj9/sczMTed/7h8Ouv2VPObv5fGIUJjjOlZrwlARPzA3cCFwDTgahGZFlZtH3Ad8Hhooaq+oqozVXUmcB7QBPwlpMo/d89X1Q0n/jbc8adN5Ryub+XW8yYhIie0jsvOHMNF00/h5y+9T8nhhhhHaIwxPYumBTAHKFHVUlVtA1YAS0IrqOoeVd0EHK8f4wrgBVVtOuFoB5hH3trLhNw0Fp7Ar/9Q/3bp6aQm+vnmM5vo7NLeFzDGmBiIJgGMAfaHvC5zyvpqKfBEWNkPRWSTiPxMRJIiLSQiN4lIsYgUV1ZWnsBm42NTWQ0b9tewfN74E/713y03I4nvXjyN9XureXr9/t4XMMaYGIgmAUTau/XpZ6qIjAKmA6tDir8FTAVmA9nANyMtq6r3qmqRqhbl5h7/FMv+9Pt3DpAY8HH5rBPJhR/1D7PGUDQuiztX76SupT0m6zTGmOOJJgGUAQUhr/OBg33czueAZ1X12J5NVcs1qBV4kGBX06DQ0dnF85vKOX9qHsOSE2KyThHhe585jaONbdz1cklM1mmMMccTTQJYB0wWkUIRSSTYlbOyj9u5mrDuH6dVgAT7Ty4DtvRxna55u7SKIw2tXDpjdEzXOz1/OFeelc+Db+ymtNIOCBtj4qvXBKCqHcCtBLtvtgNPqepWEblDRC4FEJHZIlIGXAn8RkS2di8vIuMJtiD+Frbqx0RkM7AZyAH+/eTfTv/448aDZCQF+MTUvJiv+58+PYUEv4+fvPhezNdtjDGhohoNVFVXAavCym4PmV5HsGso0rJ7iHDQWFXP60ugA0VXl/LSjgo+MTWP5AR/zNefl5HMDQsK+eXLJdzy8VpOHzM85tswxhiwK4H7bPOBWo40tHFeHH79d7tx4QSGpyTwk7/sjNs2jDHGEkAfvbLzMCL0OujbyRieksCXPz6RV3ZWUrynKm7bMcZ4myWAPnplx2HOLMgkOy0xrttZPn8cuRlJ3LnaWgHGmPiwBNAHRxpa2VhWy7lT4tf90y01McBXzp3Imt1VrN1trQBjTOxZAuiDNaXBHfGCkxz6IVpLZ49lRFoid71i1wUYY2LPEkAfrNl9lNREP9P76cyclEQ/Nyws5O/vVbKprKZftmmM8Q5LAH2wprSKs8ZlkeDvv4/t2rnjGJYc4G5rBRhjYswSQJSqGtvYWVHP2YXZ/brdjOQErps/ntVbK3ivor5ft22MGdosAUSp+0Ds3Akj+n3b159TSGqin3usFWCMiSFLAFFau7uKpICPM/Iz+33bWWmJfH7OWP64qZwDNXb/YGNMbFgCiNKG/dXMyM/s031/Y+n6BYUAPPj6ble2b4wZeiwBRKG9s4utB+s4I9+9cXnGZKZw8fRRrFi3n3q7X4AxJgYsAURh56F6Wju6mFHQ/90/ob60cAINrR08uc7uGmaMOXmWAKKw0TkHf4YL/f+hpucP5+zCbB58Yw8dnce7/bIxxvTOEkAUNu2vJSs1gYLsFLdD4UsLJ3CgpplVWw65HYoxZpCzBBCFjWU1nJGfedI3f4+F86bmMSE3jfteK0W1T7dmNsaYD7EE0Iumtg7eq6h3vf+/m88n3LCgkE1ltTZInDHmpFgC6MWOQ/V0KZw+epjboRzz2Vn5ZKclcr+dEmqMOQlRJQARWSwiO0WkRERuizB/kYi8IyIdInJF2LxOEdngPFaGlBeKyBoReV9EnnRuOD/g7CgPDr9w6qiBkwCSE/wsnV3AX7dXUFbd5HY4xphBqtcEICJ+4G7gQmAacLWITAurtg+4Dng8wiqaVXWm87g0pPxHwM9UdTJQDdxwAvHH3Y5DdaQnBcjPcv8AcKgvzB2HiPDo23vdDsUYM0hF0wKYA5SoaqmqtgErgCWhFVR1j6puAqI6N1GCR1PPA552ih4GLos66n60o7yeKadkDIgDwKFGZ6bw6dNGsmLtfprbOt0OxxgzCEWTAMYAoVcelTll0UoWkWIReVtEunfyI4AaVe3obZ0icpOzfHFlZWUfNnvyVJXth+qYekpGv243Wsvnjae2uZ3nNhxwOxRjzCAUTQKI9NO3L+cfjlXVIuDzwP+IyMS+rFNV71XVIlUtys2N343YIymvbaG+pYOpA6j/P9ScwmymnpLBQ2/usVNCjTF9Fk0CKAMKQl7nAwej3YCqHnSeS4FXgTOBI0CmiAROZJ39ZcehOgBOHaAtABHh+nPGs+NQPWvslFBjTB9FkwDWAZOds3YSgaXAyl6WAUBEskQkyZnOAc4Btmnw5+orQPcZQ8uB5/oafLxtd84A+tgATQAAS2aOITM1gYff3ON2KMaYQabXBOD0098KrAa2A0+p6lYRuUNELgUQkdkiUgZcCfxGRLY6i58KFIvIRoI7/P9S1W3OvG8C3xCREoLHBO6P5RuLhR2H6hmTmcKw5AS3Q+lRcoKfq2YXsHrrIbtXgDGmTwK9VwFVXQWsCiu7PWR6HcFunPDl3gSm97DOUoJnGA1YJYcb+NjIdLfD6NW1c8fx27+X8ru39/LNxVPdDscYM0jYlcA96OpSdh9pYELuwE8A+VmpfGraSJ5Yu4+Wdjsl1BgTHUsAPThY20xLexcTB0ECAFg+fzw1Te2s3DDgjqUbYwYoSwA9KK1sBGBCbprLkURn3oQRTBlpp4QaY6JnCaAHuyobAAZNC0BEWDZ/HNvK6yjeW+12OMaYQcASQA92VTaQkRwgJ31AjlEX0eVnjmFYcoCH7JRQY0wULAH0oLSykYm56QNuDKDjSU0M8LmiAlZvOcSh2ha3wzHGDHCWAHqwq7Jh0HT/hLp23jg6VXl8jY0Saow5PksAETS0dlBR1zpoDgCHGjcijU9MyePxtfto7bBTQo0xPbMEEEHpIDsAHG7ZvHEcaWjjhc1243hjTM8sAUSw+8jgOgU03KLJuRTmpNnBYGPMcVkCiGB/VfA2iwVZqS5HcmJ8PmHZvHFs2F/Dxv01bodjjBmgLAFEsPdoE3kZSaQk+t0O5YR99qx8UhP9PPzWHrdDMcYMUJYAIthX1cTY7MH567/bsOQEPjsrn+c3lnO0odXtcIwxA5AlgAj2D4EEALB8/jjaOrtYsW5/75WNMZ5jCSBMa0cn5XUtjB0x+BPApLwMzpk0gt+9vZeOzi63wzHGDDCWAMIcqG5GlSHRAoDgjePLa1t4cVuF26EYYwYYSwBh9jpnAA2VBHD+qSMZk5lip4QaYz4iqgQgIotFZKeIlIjIbRHmLxKRd0SkQ0SuCCmfKSJvichWEdkkIleFzHtIRHaLyAbnMTM2b+nk7B9iCcDvE66dN441u6uO3eTeGGMgigQgIn7gbuBCYBpwtYhMC6u2D7gOeDysvAlYpqqnAYuB/xGRzJD5/6yqM53HhhN8DzG172gTyQk+cjOS3A4lZq4qKiAp4OPhN218IGPMB6JpAcwBSlS1VFXbgBXAktAKqrpHVTcBXWHl76nq+870QeAwkBuTyOOk+xTQwTQKaG+y0hJZMnM0f3j3ALVN7W6HY4wZIKJJAGOA0PMIy5yyPhGROUAisCuk+IdO19DPRGRA/OQeCtcARLJ8/nia2zv53/V2SqgxJiiaBBDpp3Cf7jkoIqOAR4HrVbW7lfAtYCowG8gGvtnDsjeJSLGIFFdWVvZls32mquyvaiJ/kA4BcTynjR7O7PFZPPLWXjq77JaRxpjoEkAZUBDyOh+I+s7jIjIM+BPwHVV9u7tcVcs1qBV4kGBX00eo6r2qWqSqRbm58e09qm1up7Gtk/yslLhuxy3L5o1nX1UTf3vvsNuhGGMGgGgSwDpgsogUikgisBRYGc3KnfrPAo+o6v+GzRvlPAtwGbClL4HHw8Ga4F20RmcOzQSw+PRTGDksiYfsYLAxhigSgKp2ALcCq4HtwFOqulVE7hCRSwFEZLaIlAFXAr8Rka3O4p8DFgHXRTjd8zER2QxsBnKAf4/pOzsBB2uagaGbABL8Pq45exx/f6/y2E3vjTHeFYimkqquAlaFld0eMr2OYNdQ+HK/A37XwzrP61Ok/eBgbXcCSHY5kvhZOqeAX778Po++tZfvX3qa2+EYY1xkVwKHOFDTTKLfR07agDghKS7yMpK5ePoonl5fRkNrh9vhGGNcZAkgRHlNC6cMT8bnGzrXAESybP54Glo7+P07ZW6HYoxxkSWAEAdrmod090+3MwsyOSN/OA+/uQdVOyXUGK+yBBAimACG5gHgUCLC8nnj2VXZyBslR90OxxjjEksAjo7OLg7VtTB6+NBPAAAXnzGK7LREu2WkMR5mCcBxuL6VLh26p4CGS07wc/WcAl7aXnFsBFRjjLdYAnB8cA3A0D8G0O2as8chIvzubbswzBgvsgTgOOAkgDEeaQFAsLVzwbSRrFi3n+a2TrfDMcb0M0sAjvLa4DAQozyUACA4SmhtczvPvnvA7VCMMf3MEoDjYE0zw5IDpCdFdXH0kHF2YTbTxwznvtdK6bJRQo3xFEsAjoM1LZ45ABxKRLhxYSGlRxp5eYeNEmqMl1gCcByub2HkMO8cAA510fRRjB6ezG9fK3U7FGNMP7IE4KioayFvCN0HuC8S/D6+uKCQNbur2FRW43Y4xph+YgkA6OxSKutbPdsCALhqdgEZSQF++9put0MxxvQTSwDA0cbgRWAjh3mzBQCQkZzA0jkFrNpcfuyUWGPM0GYJADhc1wpAnodbAADXn1OIAA++bq0AY7zAEgDB/n/A011AELww7OIzRrFi3X7qWtrdDscYE2eWAIAKpwXg5S6gbl9aOIGG1g5WrN3ndijGmDiLKgGIyGIR2SkiJSJyW4T5i0TkHRHpEJErwuYtF5H3ncfykPKzRGSzs85fODeHd0VFXQsikJNuCeD0McOZN2EED76xh/bOLrfDMcbEUa8JQET8wN3AhcA04GoRmRZWbR9wHfB42LLZwPeAs4E5wPdEJMuZ/SvgJmCy81h8wu/iJB2ub2FEWhIJfmsQAXxpUSHltS08v+mg26EYY+Iomj3eHKBEVUtVtQ1YASwJraCqe1R1ExD+k/HTwIuqWqWq1cCLwGIRGQUMU9W3NHhLqkeAy072zZyoirpW6/4Jce7H8pgyMoNfvbrLhocwZgiLJgGMAfaHvC5zyqLR07JjnOle1ykiN4lIsYgUV1ZWRrnZvqmo8+5VwJH4fMIt507kvYoGXrLhIYwZsqJJAJH65qP9WdjTslGvU1XvVdUiVS3Kzc2NcrN9Yy2Aj7rkjFEUZKdw1ysldt9gY4aoaBJAGVAQ8jofiLZzuKdly5zpE1lnTHV0dnG0sZW8DGsBhAr4fdy8aCIb99fw1i67b7AxQ1E0CWAdMFlECkUkEVgKrIxy/auBC0Qkyzn4ewGwWlXLgXoRmeuc/bMMeO4E4j9pRxraULVrACK54qx8cjOSuOfVXW6HYoyJg14TgKp2ALcS3JlvB55S1a0icoeIXAogIrNFpAy4EviNiGx1lq0CfkAwiawD7nDKAG4B7gNKgF3ACzF9Z1HqvgjMqwPBHU9ygp8bFxTyeskRNu63QeKMGWqiuvuJqq4CVoWV3R4yvY4Pd+mE1nsAeCBCeTFwel+CjQe7Cvj4rpk7jrtfKeGeV0v4zbVFbodjjIkhz5/4XlHfPQ6QtQAiSU8KcN388azeWsH7FfVuh2OMiSHPJ4AjTgIYkZbociQD13XnFJKS4OdXdizAmCHF8wngaGMrWakJBOwq4B5lpyXy+bPH8tzGg+w50uh2OMaYGPH8Xu9IfZuNARSFmxdNIOATfvlyiduhGGNixBJAQ6slgCjkDUvmC3PH8ey7Zey2VoAxQ4LnE8DRxjZGpFv/fzS+/PGJJAZ8/PKl990OxRgTA55PAEfqrQUQrdyMJK6dO44/bDjArsoGt8MxxpwkTyeAlvZO6ls7yLWLwKJ288cnkhTwWyvAmCHA0wngaGMbYKeA9kVOehLL5o1j5caDlBy2VoAxg5mnE0D3NQDWBdQ3Ny2aQHKCn19YK8CYQc3bCaDBSQDWBdQnI9KTWDZvPH/cdJD37OpgYwYtTyeAow3WBXSibl40gfTEAHeu3ul2KMaYE+TpBFDptADsIHDfZaUlcvPHJ/DitgrW763qfQFjzIDj6QRwpKGV9KQAyQl+t0MZlL64oJDcjCT+64UddtcwYwYhTyeAow12EdjJSE0M8LXzJ7NuTzUv272DjRl0PJ0AbBiIk7d0dgHjRqTy4z/vpLPLWgHGDCaeTwB2APjkJPh9/L8LprCzop4/vHvA7XCMMX3g6QRwtKHNTgGNgUumj+K00cP46Yvv0drR6XY4xpgoRZUARGSxiOwUkRIRuS3C/CQRedKZv0ZExjvl14jIhpBHl4jMdOa96qyze15eLN9Ybzo6u6hqsqGgY8HnE267cCoHapp5+M09bodjjIlSrwlARPzA3cCFwDTgahGZFlbtBqBaVScBPwN+BKCqj6nqTFWdCVwL7FHVDSHLXdM9X1X79ShiVVMbqpBjB4FjYuHkXM6dkssvXyrhqHN6rTFmYIumBTAHKFHVUlVtA1YAS8LqLAEedqafBs4XEQmrczXwxMkEG0sfXARmLYBY+c7Fp9LU3slPX3zP7VCMMVGIJgGMAfaHvC5zyiLWUdUOoBYYEVbnKj6aAB50un++GyFhACAiN4lIsYgUV1ZWRhFudKq7B4KzFkDMTMrL4Nq543hi7T52HKpzOxxjTC+iSQCRdszh5/sdt46InA00qeqWkPnXqOp0YKHzuDbSxlX1XlUtUtWi3NzcKMKNTlVTMAFkpVoCiKWvnz+ZjOQEfvD8Nrs4zJgBLpoEUAYUhLzOBw72VEdEAsBwIHR8gKWE/fpX1QPOcz3wOMGupn7T3QLISkvoz80OeVlpifzfT07mjZKj/HW7XRxmzEAWTQJYB0wWkUIRSSS4M18ZVmclsNyZvgJ4WZ2ffyLiA64keOwApywgIjnOdAJwCbCFflTV2A5YCyAerpk7jom5afzwT9to6+hyOxxjTA96TQBOn/6twGpgO/CUqm4VkTtE5FKn2v3ACBEpAb4BhJ4quggoU9XSkLIkYLWIbAI2AAeA3570u+mD6qY2MpIDJPg9fSlEXCT4fXznkmnsOdrE/a/vdjscY0wPAtFUUtVVwKqwsttDplsI/sqPtOyrwNywskbgrD7GGlNVjW1k21XAcfOJKXl8atpIfvHS+1w6czRjMlPcDskYE8azP3+rm9qs+yfOvveZaSjKv63c6nYoxpgIPJ0ArAUQX/lZqXzt/Mn8ZVsFL22vcDscY0wY7yaAxnZrAfSDGxdMYGJuGt9buZXmNhsnyJiBxLMJIHgMwE4BjbfEgI8fXHY6ZdXN3PNqidvhGGNCeDIBNLd10tzeSZZ1AfWL+RNzuPzMMfz6b7soOWw3kTdmoPBkAqh2rgLOti6gfvPti04lLSnAvzy9yW4cY8wA4ckEUOVcBZxpCaDf5GYkcfsl03hnX40NGW3MAOHJBHCsBWBdQP3q8jPH8Ikpudy5eif7jja5HY4xnufJBNDdArCDwP1LRPjh5dPx+4Tbfr/JBoszxmWeTADHBoKzLqB+NzozhW9fdCpv7jrKinX7e1/AGBM33kwATe2IwPAUawG44eo5BcybMIIf/mk7ZdXWFWSMWzyaANoYnpJAwAaCc4WI8OMrzgDgG09ttLOCjHGJJ/eAVY1tdgqoywqyU/n+paexdncV9/69tPcFjDEx58kEUN3UZheBDQCfnTWGi6afwk9f3MmWA7Vuh2OM53gyAVTZOEADgojww8umk52WyNdXvGtjBRnTzzyZAKptHKABIystkf++cga7Khv5zxe2ux2OMZ7iuQSgqlTZvQAGlIWTc7lhQSGPvLWXFzaXux2OMZ7huQTQ0t5FW0eXDQMxwHxz8VRmFGTyL09vYu/RRrfDMcYTokoAIrJYRHaKSImI3BZhfpKIPOnMXyMi453y8SLSLCIbnMevQ5Y5S0Q2O8v8QkQkVm/qeGqagxeB2TUAA0tiwMddV5+Jzyd85bF3aGm34wHGxFuvCUBE/MDdwIXANOBqEZkWVu0GoFpVJwE/A34UMm+Xqs50Hl8OKf8VcBMw2XksPvG3Eb3a5nYAMlMtAQw0Bdmp/PRzM9h6sI4fPL/N7XCMGfKiaQHMAUpUtVRV24AVwJKwOkuAh53pp4Hzj/eLXkRGAcNU9S0NDgjzCHBZn6M/ATVNwQRgLYCB6fxTR3Lzxyfw2Jp9PLfhgNvhGDOkRZMAxgChg7aUOWUR66hqB1ALjHDmFYrIuyLyNxFZGFK/rJd1AiAiN4lIsYgUV1ZWRhHu8XW3ACwBDFz/dMEU5ozP5pvPbLLrA4yJo2gSQKRf8uHX7vdUpxwYq6pnAt8AHheRYVGuM1ioeq+qFqlqUW5ubhThHl+ttQAGvAS/j7uvmUV2aiI3PVLMkYZWt0MyZkiKJgGUAQUhr/OBgz3VEZEAMByoUtVWVT0KoKrrgV3Ax5z6+b2sMy7sGMDgkJuRxL3LiqhqauOW362nraPL7ZCMGXKiSQDrgMkiUigiicBSYGVYnZXAcmf6CuBlVVURyXUOIiMiEwge7C1V1XKgXkTmOscKlgHPxeD99KqmuQ2/T0hPCvTH5sxJOH3McO68Ygbr9lTzvZVb7f4BxsRYr3tBVe0QkVuB1YAfeEBVt4rIHUCxqq4E7gceFZESoIpgkgBYBNwhIh1AJ/BlVa1y5t0CPASkAC84j7irbW5neEoC/XTWqTlJn5kxmu3lddzz6i4m56XzxQWFbodkzJAR1c9gVV0FrAoruz1kugW4MsJyzwDP9LDOYuD0vgQbCzVN7db/P8j80wVT2FXZwA/+tI3RmcksPn2U2yEZMyR47krg7haAGTx8PuHnS8/kzIJMvr5iA8V7qnpfyBjTK0sAZlBITvBz3/LZjM5M4cZHitlV2eB2SMYMep5MAHYG0OCUnZbIw9fPIeATlj+wloq6FrdDMmZQ81wCsGMAg9vYEancv3w21Y1tXHPfGo7aNQLGnDBPJYCuLqWupZ1MSwCD2oyCTO6/bjb7q5q49v61xy7uM8b0jacSQH1rB6owzBLAoDd3wgh+c+1ZvH+4nuUPrqWhtcPtkIwZdDyVALp/Kdq9AIaGc6fkcdfnZ7H5QC1ffGgdjZYEjOkTbyUAGwhuyPn0aafws6tmsn5vNcsfWEtdi3UHGRMtTyWA7pvB2FlAQ8ulM0Zz19VnsmF/DV+4bw01TW1uh2TMoOCpBGAtgKHrwumjuHfZWew4VM/Se9+2EUSNiYKnEkD3zWDsLKCh6bypI3lg+Wz2Hm3ic79+i/1VTW6HZMyA5qkE0N0CsLOAhq4Fk3N49IY5HG1s4/J73mRzmd1QxpieeC4BJAV8JCf43Q7FxFHR+GyeuWU+SQEfV937Fq/sPOx2SMYMSN5KAE02DIRXTMpL59mvzKcwJ40bHy7mibX73A7JmAHHUwmgprnNDgB7SN6wZJ68eR4LJuXwrd9v5vbnttDeaXcWM6abpxJAbXM7mSl2EZiXpCcFuH95ETctmsAjb+3lCzZ+kDHHeCoB1DS12wFgDwr4fXz7olP5n6tmsmF/DZfe9QZbDtjBYWM8lQDqWzqsC8jDLjtzDM/cMh9V5R/ueZOH3tht9xk2nhZVAhCRxSKyU0RKROS2CPOTRORJZ/4aERnvlH9KRNaLyGbn+byQZV511rnBeeTF6k31pK6lnYxkuxm8l50+ZjjPf20hCybn8P0/buPmR9fblcPGs3pNACLiB+4GLgSmAVeLyLSwajcA1ao6CfgZ8COn/AjwGVWdDiwHHg1b7hpVnek84nquXleX0tDawTBLAJ6XnZbI/cuL+M7Fp/LKzsNc9PPXWGe3mTQeFE0LYA5QoqqlqtoGrACWhNVZAjzsTD8NnC8ioqrvqupBp3wrkCwiSbEIvK8a2mwoaPMBEeHGhRN45pb5BPw+Pvebt/j357fR0t7pdmjG9JtoEsAYYH/I6zKnLGIdVe0AaoERYXU+C7yrqqGnYDzodP98V0Qk0sZF5CYRKRaR4srKyijCjay+JThUsHUBmVBn5Gey6usL+QT8ykIAAA6QSURBVPycsdz3+m4u+vlrrN9b7XZYxvSLaBJApB1z+JGz49YRkdMIdgvdHDL/GqdraKHzuDbSxlX1XlUtUtWi3NzcKMKNrK57GIhkawGYD0tPCvDDy6fzuxvOprWjiyt//SY//NM2u7+AGfKiSQBlQEHI63zgYE91RCQADAeqnNf5wLPAMlXd1b2Aqh5wnuuBxwl2NcXNBy0ASwAmsgWTc/jzPy7kqtlj+e1ru/nkT//GC5vL7UwhM2RFkwDWAZNFpFBEEoGlwMqwOisJHuQFuAJ4WVVVRDKBPwHfUtU3uiuLSEBEcpzpBOASYMvJvZXjO9YCSLEuINOzjOQE/vMfpvPMLfMYnpLALY+9w/IH17HnSKPboRkTc70mAKdP/1ZgNbAdeEpVt4rIHSJyqVPtfmCEiJQA3wC6TxW9FZgEfDfsdM8kYLWIbAI2AAeA38byjYWrbw0mAGsBmGicNS6b5//PAm6/ZBrv7K3mgp/9nX9/fpudMmqGFBlMzduioiItLi4+oWUffnMP31u5leLvfJKcdFdORDKDVEVdC/+9eidPv1NGRlKAW8+bxLJ5421UWTNoiMh6VS0KL/fMlcD1Ld0tAOsCMn0zclgyd145gxe+vpBZ47L4j1U7OP8nf+OpdfttcDkzqHkoAXSQFPCRFLBfbebETD1lGA9dP4fHbjybEemJ/Mszmzj3zld5bM1eWjvs+gEz+HgmAdS12EBwJjbOmZTDc189hwevm03esCT+9dktnHvnqzz4xm47ddQMKh5KAB3W/WNiRkT4xNQ8fn/LfH53w9kUZKXyb3/cxtz/fIn/WLWdsmq7H7EZ+DyzR6xrbreLwEzMiQgLJuewYHIO7+yr5oHXd3P/67u577VSFp9+Csvmjefswmx6uNDdGFd5JgHUWwvAxNmssVnM+nwWB2uaeeStvTyxdh+rNh+iMCeNK4vyuWJWPnnDkt0O05hjPNQFZMcATP8YnZnCbRdO5e1vnc9PrpxBbnoSP/7zTub918vc+HAxf95SboPOmQHBMz+J61tsKGjTv1IS/Xz2rHw+e1Y+uyobeKp4P8+sP8Bft1eQnhTgU9NG8pkZo1gwKZfEgGd+i5kBxDN7xPqWdrsK2LhmYm4637rwVP75gim8VXqU5zeW88KWcp599wDDUxL41LSRfPLUPBZMziU9yTP/lsZlnvimtXV00dLeZS0A47qA38fCybksnJzLDy47nddLKvnjxnL+svUQT68vI9Hv4+wJ2Zw/NY/zTx1JQXaq2yGbIcwTe8QPrgK2FoAZOBIDPs6bOpLzpo6kvbOL9XureWl7BS/tOMz3/7iN7/9xG+NGpDJ/4gjmTcxh3oQR5GbYMCYmdjyRAOqcoaBtJFAzUCX4fcydMIK5E0bwrxdPY/eRRl7deZg3SoLdRU+sDd6T6WMj05k7YUTwjKOxWRRkp9gppuaEeWKPeKwFkGQtADM4FOakUZhTyPXnFNLR2cXWg3W8uesob+46wtPry3jkrb0A5KQnMrMgi1njMplZkMlpo4cz3M52M1HyRAKoa+5uAdg/hhl8An4fMwoymVGQyS3nTqSjs4udFfW8u6/GeVTz1+0Vx+rnZ6UwbdQwpo0exqmjhjFt1DDys6ylYD7KEwnARgI1Q0nA7+O00cM5bfRwvjB3HADVjW1sLKthW3kd28vr2Xawlhe3V9A92nt6UoAJuWlMzE1nQk4aE/PSmZCbxvgRaTastYd5Yo9oN4Q3Q11WWiLnTsnj3Cl5x8qa2jrYeaie7eX17DxUR+mRRtburuLZdw8cqyMSbDEUZKWSn5VCflYqBdnB5/ysFEZmJOPzWcthqPLEHrGupft2kNYFZLwjNTHAmWOzOHNs1ofKm9o6KK1spPRII7sON1B6pJGy6iZe2VlJZX3rh+om+IXRmcFEkDcsiZHDksnLcJ6HJZGXkczIYUmkJwWsi2kQ8kgC6EAE0hM98XaNOa7UxACnjxnO6WOGf2ReS3snB2qa2V/VRFl1s/No4nB9K1sO1PLS9sM0RxjGIjnBR3ZqIllpiWSnJZKV+sFzVlrCsdfDUxIYlpxARnKA9OQACX67AtpNUe0RRWQx8HPAD9ynqv8VNj8JeAQ4CzgKXKWqe5x53wJuADqBr6nq6mjWGUt1ze2kJwWsKWtML5IT/EzMTWdibnrE+apKQ2sHFXWtHK5v4XBdKxV1LVTWt1Ld1E51UxtVjW3sr2qiqrHt2CnYPW/PR4aTEDKSE8hICjjTAdKTEkhL8pOc4Cc10U9Kgp8U5zk1MUBKoo+UhMCxsu5nG1Yjer0mABHxA3cDnwLKgHUislJVt4VUuwGoVtVJIrIU+BFwlYhMA5YCpwGjgb+KyMecZXpbZ8wExwGy7h9jTpaIODvsBCblRU4Sodo7u6hxEkN1YxvVTe00tHZQ39JOfUvIc2vHsdeH6lpocKab2jvp623L/T4hwS8k+n0kBvwkBXwkBnzBsoDPKQ/OSwwrSwh59vuEgE9Cnn0fvPZHLk/wh9ULqe+TkIePY9N+X/Bz9YngF0EEfL7gtE+C8/w+YXhKAv4Y/4iNpgUwByhR1VIAEVkBLAFCd9ZLgO87008Dd0mwQ3AJsEJVW4HdIlLirI8o1hkzwXGArPvHmP6W4PeRm5F0wlcwqyqtHV00t3XS1N5Jc1snLe2dNLV10tzeSXNbh/PcRVNbBy3twfK2ji7aO4PLtnV00dbZRbvz3OaU1Ta3f6Ss3Zlu7+qis0tp7+xj9omjv37j41El3b6IZq84Btgf8roMOLunOqraISK1wAin/O2wZcc4072tEwARuQm4CWDs2LFRhPtRMwoymRjjD84YE38iQnJCsBsoq/fqcdHVpXR0KZ1dSoeTGD54rXR29lDe1UVHZ/B1e5fSpYqq0tnFR6aPPT702pnuCk7npsd+GJBoEkCkNkd4WuypTk/lkTrpIqZaVb0XuBegqKjohNLxVz8x6UQWM8YYfD4h8VjXy9C6ZiKaoyVlQEHI63zgYE91RCQADAeqjrNsNOs0xhgTR9EkgHXAZBEpFJFEggd1V4bVWQksd6avAF5WVXXKl4pIkogUApOBtVGu0xhjTBz12gXk9OnfCqwm2P55QFW3isgdQLGqrgTuBx51DvJWEdyh49R7iuDB3Q7gq6raCRBpnbF/e8YYY3oi2tdzrFxUVFSkxcXFbodhjDGDioisV9Wi8HK7YsIYYzzKEoAxxniUJQBjjPEoSwDGGONRg+ogsIhUAntPcPEc4EgMw4kVi6tvLK6+sbj6ZqjGNU5Vc8MLB1UCOBkiUhzpKLjbLK6+sbj6xuLqG6/FZV1AxhjjUZYAjDHGo7yUAO51O4AeWFx9Y3H1jcXVN56KyzPHAIwxxnyYl1oAxhhjQlgCMMYYjxpyCUBEFovIThEpEZHbIsxPEpEnnflrRGR8P8RUICKviMh2EdkqIl+PUOdcEakVkQ3O4/Z4x+Vsd4+IbHa2+ZGR9iToF87ntUlEZvVDTFNCPocNIlInIv8YVqdfPi8ReUBEDovIlpCybBF5UUTed54j3qxKRJY7dd4XkeWR6sQ4rjtFZIfzd3pWRDJ7WPa4f/M4xPV9ETkQ8re6qIdlj/u/G4e4ngyJaY+IbOhh2Xh+XhH3Df32HVPn1mRD4UFwaOldwAQgEdgITAur8xXg1870UuDJfohrFDDLmc4A3osQ17nA8y58ZnuAnOPMvwh4geDd3eYCa1z4mx4ieCFLv39ewCJgFrAlpOzHwG3O9G3AjyIslw2UOs9ZznRWnOO6AAg40z+KFFc0f/M4xPV94J+i+Dsf93831nGFzf8JcLsLn1fEfUN/fceGWgvg2A3sVbUN6L7ZfKglwMPO9NPA+SIS6daVMaOq5ar6jjNdD2zng3sjD3RLgEc06G0gU0RG9eP2zwd2qeqJXgF+UlT17wTvcREq9Dv0MHBZhEU/DbyoqlWqWg28CCyOZ1yq+hdV7XBevk3wTnv9qofPKxrR/O/GJS7n//9zwBOx2l60jrNv6Jfv2FBLAJFuYB++o/3QDeyB7hvY9wuny+lMYE2E2fNEZKOIvCAip/VTSAr8RUTWi8hNEeZH85nG01J6/sd04/MCGKmq5RD8BwbyItRx+3P7IsGWWyS9/c3j4Vana+qBHroz3Py8FgIVqvp+D/P75fMK2zf0y3dsqCWAk7mBfdyJSDrwDPCPqloXNvsdgt0cM4BfAn/oj5iAc1R1FnAh8FURWRQ2383PKxG4FPjfCLPd+ryi5ebn9q8E78D3WA9Vevubx9qvgInATKCcYHdLONc+L+Bqjv/rP+6fVy/7hh4Xi1DWp89sqCWAk7mBfVyJSALBP/Bjqvr78PmqWqeqDc70KiBBRHLiHZeqHnSeDwPPEmyKh4rmM42XC4F3VLUifIZbn5ejorsbzHk+HKGOK5+bcyDwEuAadTqKw0XxN48pVa1Q1U5V7QJ+28P23Pq8AsA/AE/2VCfen1cP+4Z++Y4NtQRwMjewjxunj/F+YLuq/rSHOqd0H4sQkTkE/zZH4xxXmohkdE8TPIi4JazaSmCZBM0Farubpv2gx19mbnxeIUK/Q8uB5yLUWQ1cICJZTpfHBU5Z3IjIYuCbwKWq2tRDnWj+5rGOK/SY0eU9bC+a/914+CSwQ1XLIs2M9+d1nH1D/3zH4nFk280HwbNW3iN4RsG/OmV3EPynAEgm2KVQAqwFJvRDTAsINs02ARucx0XAl4EvO3VuBbYSPPvhbWB+P8Q1wdneRmfb3Z9XaFwC3O18npuBon76O6YS3KEPDynr98+LYAIqB9oJ/uK6geAxo5eA953nbKduEXBfyLJfdL5nJcD1/RBXCcE+4e7vWPfZbqOBVcf7m8c5rked784mgju2UeFxOa8/8r8bz7ic8oe6v1Mhdfvz8+pp39Av3zEbCsIYYzxqqHUBGWOMiZIlAGOM8ShLAMYY41GWAIwxxqMsARhjjEdZAjDGGI+yBGCMMR71/wH2oCWCIhi/jAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.arange(0, 20, 0.001)\n",
    "y = f_gamma(x)\n",
    "plt.plot(x, y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Функция распределения"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAdQ0lEQVR4nO3deXRddb338fc3c9ukY9JQkoa0tHSgQAuxZagKKAJVKCggBQV8wPo8ylWXXr0oV/Q6LJd6vXp9Lo9aBAEtQwGBqmVwQAZpCymldKKQhjZDh6RtxmY+5/v8kdMS2qQ5bU+yT875vNY66+zhd06+a3P6WZvf3vv3M3dHRESGvpSgCxARkdhQoIuIJAgFuohIglCgi4gkCAW6iEiCSAvqD+fm5npxcXFQf15EZEhas2bNHnfP621fYIFeXFxMaWlpUH9eRGRIMrPtfe1Tl4uISIJQoIuIJAgFuohIglCgi4gkCAW6iEiC6DfQzeweM6sxsw197Dcz+4WZlZnZG2Z2ZuzLFBGR/kRzhn4vcMkR9l8KTI28FgO/PP6yRETkaPV7H7q7v2BmxUdoshC437vH4V1lZqPNbIK774xRjSKSZMJhpyMUJhR2usJO+MC7O6Hwu69otx18efd3hR3C7jjg7viBdQene5kebQ7uO7AePvDZSNsebcIOzrvtw/5uO6f7Qx+akc8ZE0fH/LjF4sGiAqCyx3pVZNthgW5mi+k+i6eoqCgGf1pEBko47LR0hmhq62R/e4i2zu5Xa2eI1o4QbV1h2jq61w9u7wzR1hGirTNMW1eIzlCYjq4wHSGnsytMZ6j71RFyOrpCdIb83W1dYTpCYTpD3eGbqMwgf1RW3Aa69bKt1/8a7r4EWAJQUlKSuP/FROJEW2eIupYO9u3voL6lM/Lewb79nTS2ddLc1kVzexdN7V00t3XSFFlvbuuiuaOLo53/JjMthWEZqQxLTyUzLYWMtBTSU7tfGakpZKankJ2VRkZqCulp3dvSU+3dNge3pZCWaqSnGilmpKUYqSlGSkr3cooZaZF9qX1sS00xUnvZdmDZADPDrDvEUiLLKdYdaSmRNge2d7czUizyucg+jIPbUiJtDrRPsUO/o7e4jJ1YBHoVMLHHeiGwIwbfKyK96AyF2dPczu7GdnY3tlHT1E5NY9vB5dqmdur2d1DX0klrZ6jP7xmRkUp2VhrZmWlkZ6UzMiuN/JFZkfU0cg68Z6UzPBLSWempBwP7wHJWJMSz0lJJSRnYwJIji0WgLwduNbOHgHlAg/rPRY5dR1eYHfWtVNa1UFXXSuW+7vequhYq61rZ09x+2JlzikFudib5I7MYn5PJ9BNGMmZ4OmNGZDB2REb38vAMxozIYMzwDEYPTyc9VXctJ5p+A93MHgTOB3LNrAr4NpAO4O6/AlYAC4AyoAX4zEAVK5JI9u3vYGttM1trmimrae5ert1PZV3LewI7LcU4cfQwCscM48Jp45kwOovxOVnkj8w8+D4uO5NUnR0nvWjuclnUz34HvhCzikQSTCjslNc2s3FHI5t2NrJxRwObdzaxb3/HwTaZaSlMzsvmjImjuWJOAUVjh1M4ZhgTxw4nPyeTNJ1NSxQCGz5XJFHtqG+ldHsdr22v4/XKet7c1UhbZxiAjLQUpp+Qw0dm5jNlfDZTxmdzcl42BaOHqf9ZjpsCXeQ4uDtbdjexcuvegyG+s6ENgGHpqZxWOIrr5p7EqSeO5NSCkZycl62+axkwCnSRo7SroY2Xyvbw0tu1vFS2lz3N7QAUjB5GSfFYzioazVknjWXGhBx1lcigUqCL9MPdWV/dwLMbd/OXTbvZsrsJgHEjMjhvSi7zp+Zy3pRcCkYPC7hSSXYKdJFehMLO6vK9PL1xF3/ZtJudDW2kGMydNJZvLpjO/Cl5TD8hR/3eElcU6CIR7s7GHY08sbaa5et2UNPUTlZ6Ch88JY+vfmQaF04fz9gRGUGXKdInBbokvb3N7TyyporH1lTxdk0z6anG+dPGc8XsAi6cPp5hGalBlygSFQW6JCV359VtdSxdvZ2n1u+iIxSm5KQx/ODKWSyYNYExOhOXIUiBLkmlvSvEk2t3cPdL77BldxM5WWlcN6+I6+cVMTU/J+jyRI6LAl2SQmNbJw+sruCel96hpqmdGRNG8qNPnMZlZ5zI8Az9M5DEoF+yJLSGlk6WvLiV+1/eTlN7F/On5PLTa85g/pTcAR/KVGSwKdAlITW3d/Hbl95hyYvlNLV18dHTJvC/P3gypxWOCro0kQGjQJeE0tEV5nertnPnc2Xs29/BRTPz+cpFpzBjwsigSxMZcAp0SRjPbanhe3/aRHntfuZPyeVfL57G7AGY5kskXinQZcgrr23me3/axHNbapmcO4Lf3vQ+Lpg+PuiyRAadAl2GrI6uML9+fiv/9+9lZKSlcPuCGdx4bjEZaRoQS5KTAl2GpHWV9fzbY2/w5q4mPnr6BL592UzG52QFXZZIoBToMqS0d4X4r2ff4q4Xy8nLyeSuG0q4aGZ+0GWJxAUFugwZZTVNfPHB19m0s5FFc4v4xoLpjMxKD7oskbihQJe45+4sXV3B9/+8ieEZadx9YwkfmqGzcpFDKdAlrjW1dfK1R97g6Y27+MApefzn1aerr1ykDwp0iVtv727ic79fw/a9Ldy+YAY3z5+kCSVEjkCBLnHpqfU7+ddH1jEsI5Wlt8zj7Mnjgi5JJO4p0CWuuDs/+8tb/OLvZcwpGs0vrz+LE0api0UkGgp0iRvtXSFue2w9j6+t5pqSQr53xSwy0zRbkEi0FOgSF+pbOlj8uzW88s4+vnbxND5//ska3lbkKCnQJXDV9a3ccPdqKve18t/Xzmbh7IKgSxIZkhToEqhte/Zz/W9W09jayf03z9XFT5HjoECXwGzZ1cSn7l5NKOw8uPhsZhVo8gmR46Fh6SQQ66sa+OSSlaQYPKwwF4kJnaHLoNtQ3cD1v1nFyGHpPHDL2RSNGx50SSIJIaozdDO7xMy2mFmZmd3Wy/4iM3vOzNaa2RtmtiD2pUoieGt3E5++ezU5Wek8tFhhLhJL/Qa6maUCdwKXAjOBRWY285Bm/w4sc/c5wLXA/4t1oTL0ldc2c91dq8lIS+GBz86jcIzCXCSWojlDnwuUuXu5u3cADwELD2njwIFZeEcBO2JXoiSCyn0tXP+b1d0jJ95yNieNGxF0SSIJJ5pALwAqe6xXRbb19B3gU2ZWBawA/qW3LzKzxWZWamaltbW1x1CuDEV1+zu48bevsL+9i9/fMo8p47ODLkkkIUUT6L09rueHrC8C7nX3QmAB8DszO+y73X2Ju5e4e0leXt7RVytDTltniFvuL6WqrpXf3Pg+ZkwY2f+HROSYRBPoVcDEHuuFHN6lcjOwDMDdVwJZQG4sCpShKxR2vvzQ67xWUcfPPzmbuZPGBl2SSEKLJtBfBaaa2SQzy6D7oufyQ9pUAB8CMLMZdAe6+lSS3Pf+tImnN+7iWx+dyYLTJgRdjkjC6zfQ3b0LuBV4BthM990sG83su2Z2eaTZV4HPmtk64EHgJnc/tFtGksgDqyu49+Vt3Dx/Ev9r/qSgyxFJClE9WOTuK+i+2Nlz2x09ljcB58W2NBmqXnlnH3c8uYEPnpLHNxfMCLockaShR/8lpqrrW/k/v19D0djh/GLRHFI1ZZzIoFGgS8y0doRYfH8pHV1hltxQwqhh6UGXJJJUNJaLxIS7c/sT69m0s5G7byzRveYiAdAZusTEI2uq+MNr1XzxwqlcOD0/6HJEkpICXY7bll1N3PHkBs49eRxf/NDUoMsRSVoKdDku+9u7+PzSNWRnpvPza2frIqhIgNSHLsflW09soHzPfpbePI/xOVlBlyOS1HSGLsfsyder+cPa7n7zc6dopAeRoCnQ5ZjsbGjlW09sYE7RaP7lwilBlyMiKNDlGITDztcffYPOkPOza2aTlqqfkUg80L9EOWr3r9zGi2/v4d8/NoPiXE1UIRIvFOhyVMpqmvjhU29ywbQ8rptbFHQ5ItKDAl2iFgo7X3v0DYZnpPKjq07HTLcoisQTBbpE7f6V21hbUc+3LztVtyiKxCEFukSlqq6FnzyzhfOn5bFw9olBlyMivVCgS7/cndsf3wDA96+Ypa4WkTilQJd+Pfn6Dp5/q5avXzyNwjHDgy5HRPqgQJcjqtvfwX/8cSNnFo3m0+cUB12OiByBAl2O6MfPbKGxrYsffvx0DbwlEucU6NKnN6rqeejVCm46t5hpJ+QEXY6I9EOBLr0Kh51vPbmR3OxMvvxhjXEuMhQo0KVXy0orWVdZzzcXTCcnS3ODigwFCnQ5TH1LBz96+k3mFo/litkFQZcjIlFSoMthfvrsWzS2dfEfC0/VPeciQ4gCXd6jrKaJB16p4Pp5RcyYMDLockTkKCjQ5T1+uOJNhqen8iVN9iwy5CjQ5aCXy/bwtzdr+MKFUxiXnRl0OSJylBToAnTfpviDFZspGD2Mm84tDrocETkGCnQB4PG11Wzc0cjXL5lGVnpq0OWIyDFQoAutHSH+89ktnF44istO19C4IkNVVIFuZpeY2RYzKzOz2/poc42ZbTKzjWb2QGzLlIH025ffYWdDG7cvmEGKxmsRGbLS+mtgZqnAncBFQBXwqpktd/dNPdpMBb4BnOfudWY2fqAKlthqaO3k18+Xc8G0POZNHhd0OSJyHKI5Q58LlLl7ubt3AA8BCw9p81ngTnevA3D3mtiWKQPl7pfeoaG1k69+ZFrQpYjIcYom0AuAyh7rVZFtPZ0CnGJm/zSzVWZ2SW9fZGaLzazUzEpra2uPrWKJmX37O7j7xXIunXUCswpGBV2OiBynaAK9t05VP2Q9DZgKnA8sAn5jZqMP+5D7EncvcfeSvLy8o61VYuxXz2+lpTPEVy46JehSRCQGogn0KmBij/VCYEcvbZ509053fwfYQnfAS5yqaWzjvpe3ceXsAqbma6xzkUQQTaC/Ckw1s0lmlgFcCyw/pM0TwAUAZpZLdxdMeSwLldj6n+fKCIWdL2msc5GE0W+gu3sXcCvwDLAZWObuG83su2Z2eaTZM8BeM9sEPAd8zd33DlTRcnyq61t58JUKri6ZyEnjRgRdjojESL+3LQK4+wpgxSHb7uix7MBXIi+Jc7/6x1YAbr1wSsCViEgs6UnRJLO7sY2HSyv5xJmFFIweFnQ5IhJDCvQkc9cL5YTCzufP19m5SKJRoCeRvc3tLF1dwcIzTqRo3PCgyxGRGFOgJ5F7/vkObV0hPn/ByUGXIiIDQIGeJBpaOrnv5e0smDWBKeN137lIIlKgJ4n7Vm6jub1Ld7aIJDAFehLY397FPf98hw/PyNfEzyIJTIGeBJaVVlLf0qm+c5EEp0BPcF2hMHe/9A7vKx7DmUVjgi5HRAaQAj3BPbVhF1V1rXz2/ZODLkVEBpgCPYG5O0teKGdy7gg+PCM/6HJEZIAp0BPYyvK9rK9u4Jb3T9ZcoSJJQIGewO56oZzc7Aw+fuahE0yJSCJSoCeot3Y38dyWWm44p5is9NSgyxGRQaBAT1B3vVBOVnoKnz77pKBLEZFBokBPQDWNbTzxejXXlExkzIiMoMsRkUGiQE9AS1dX0BV2PnPepKBLEZFBpEBPMO1dIZauruCCaeOZlKvp5USSiQI9waxYv5M9ze3ceG5x0KWIyCBToCeYe1/ezuS8Ebx/Sm7QpYjIIFOgJ5C1FXWsq6znpnOL9SCRSBJSoCeQ+17eRk5mGh8/szDoUkQkAAr0BFHT2Maf1+/kqpJCsjPTgi5HRAKgQE8QB25VvPGc4qBLEZGAKNATQEdX+OCtisW6VVEkaSnQE4BuVRQRUKAnhN+v2s6kXN2qKJLsFOhD3Ju7GindXsf184p0q6JIklOgD3EPrK4gIy2FT+hWRZGkp0Afwva3d/GH16r52GkTNKqiiCjQh7I/rttBc3sX180rCroUEYkDUQW6mV1iZlvMrMzMbjtCu6vMzM2sJHYlSl+Wrq5gWn4OZ500JuhSRCQO9BvoZpYK3AlcCswEFpnZzF7a5QBfBFbHukg53BtV9ayvbuD6s4sw08VQEYnuDH0uUObu5e7eATwELOyl3feAHwNtMaxP+vDA6gqGpadyxRxNAC0i3aIJ9AKgssd6VWTbQWY2B5jo7n860heZ2WIzKzWz0tra2qMuVro1tnXy5Os7WDj7REZmpQddjojEiWgCvbf/n/eDO81SgJ8BX+3vi9x9ibuXuHtJXl5e9FXKezyxtprWzhDXz9ME0CLyrmgCvQqY2GO9ENjRYz0HmAX8w8y2AWcDy3VhdGC4O0tXVXB64ShOKxwVdDkiEkeiCfRXgalmNsnMMoBrgeUHdrp7g7vnunuxuxcDq4DL3b10QCpOcmu217FldxPXzdWtiiLyXv0Gurt3AbcCzwCbgWXuvtHMvmtmlw90gfJeD7xSQXZmGpedcWLQpYhInIlqJgR3XwGsOGTbHX20Pf/4y5LeNLV1smL9Tq6cU8gITWIhIofQk6JDyB/X7aStM8wn3zex/8YiknQU6EPIw6WVTMvP4QxdDBWRXijQh4gtu5pYV1nPNe+bqCdDRaRXCvQhYllpJempxpV6MlRE+qBAHwI6usI8vraai2bmM1bD5IpIHxToQ8BfN+9m3/4OrinRxVAR6ZsCfQh4+NVKJozK4v1TNVyCiPRNgR7ndtS38sLbtVx1ViGpmjNURI5AgR7nHltThTtcfZa6W0TkyBTocSwcdpatqeTck8dRNG540OWISJxToMexVeV7qdzXqidDRSQqCvQ49nBpJTlZaVx86glBlyIiQ4ACPU41tHTy1IZdXDG7gKz01KDLEZEhQIEep5avq6ajSwNxiUj0FOhx6uHSSmZOGMmsAg3EJSLRUaDHoY07GthQ3aizcxE5Kgr0OPRIaRUZaSksnK1ZiUQkegr0ONPWGeLxtdVcfOoJjB6ugbhEJHoK9Djz7KbdNLR28kkNxCUiR0mBHmceKa2kYPQwzj15XNCliMgQo0CPI1V1LbxUtoerSwpJ0UBcInKUFOhx5NE1VQBcdVZhwJWIyFCkQI8T4bDzSGkV86fkUjhGA3GJyNFToMeJleV7qa5v5WpdDBWRY6RAjxPLSisZmZXGR2bmB12KiAxRCvQ4cHAgrjkaiEtEjp0CPQ4cGIhLk0CLyPFQoMeBZaVVGohLRI6bAj1gm3Y0sr66gWtKdKuiiBwfBXrAHllTSUZqCgtnFwRdiogMcVEFupldYmZbzKzMzG7rZf9XzGyTmb1hZn8zs5NiX2riae8K8cTaai46NZ8xIzQQl4gcn34D3cxSgTuBS4GZwCIzm3lIs7VAibufDjwK/DjWhSaiv26qoa5FA3GJSGxEc4Y+Fyhz93J37wAeAhb2bODuz7l7S2R1FaAO4SgsK63kxFFZnDclN+hSRCQBRBPoBUBlj/WqyLa+3Aw81dsOM1tsZqVmVlpbWxt9lQloR30rL7xdy1VnFZKqgbhEJAaiCfTe0sZ7bWj2KaAE+Elv+919ibuXuHtJXl5e9FUmoD+8VoU7XHWWultEJDbSomhTBfRMnUJgx6GNzOzDwO3AB929PTblJaZw2FlWWsU5k8dRNE4DcYlIbERzhv4qMNXMJplZBnAtsLxnAzObA/wauNzda2JfZmJZVb6Xin0tmgRaRGKq30B39y7gVuAZYDOwzN03mtl3zezySLOfANnAI2b2upkt7+PrBFj6SgWjh6dzyawTgi5FRBJINF0uuPsKYMUh2+7osfzhGNeVsPY0t/Psxl3ccE6xBuISkZjSk6KD7NE1VXSGnEVz1d0iIrGlQB9E4bDz4CsVzJ00linjc4IuR0QSjAJ9EK0s38v2vS1cN7co6FJEJAEp0AfRA6t1MVREBo4CfZDUNrXzzMZdfOLMQl0MFZEBoUAfJI+uqaIr7CxSd4uIDBAF+iB478XQ7KDLEZEEpUAfBP94q4aKfS18+mwNEy8iA0eBPgh++89t5I/M1MVQERlQCvQBVlbTzItv7+FT804iPVWHW0QGjhJmgN2/chsZqSksmqeLoSIysBToA6iprZPH1lTxsdMnkJudGXQ5IpLgFOgD6NE1VezvCHHjucVBlyIiSUCBPkBCYee+l7cxp2g0Z0wcHXQ5IpIEFOgD5JmNu9i2t4Vb5k8OuhQRSRIK9AHg7vzq+a0UjxuuWxVFZNAo0AfAyq17eaOqgc9+YDKpKb3NsS0iEnsK9AHwy+e3kpudySfOLAy6FBFJIgr0GNtQ3cCLb+/hM+dpijkRGVwK9Bj7+V/fJicrjU9p3BYRGWQK9BhaV1nPXzfvZvH7JzNqWHrQ5YhIklGgx9BP//IWY4an85n5k4IuRUSSkAI9Rl7dto8X3qrlcx88mezMtKDLEZEkpECPgXDY+f6fNzM+J5MbzlHfuYgEQ4EeA4+vrWZdZT3/dsl0hmfo7FxEgqFAP07N7V386Ok3OWPiaK6cUxB0OSKSxBTox+nHT79JbXM737lsJil6KlREAqRAPw4rt+7l/pXbuencYuYUjQm6HBFJcgr0Y1Tf0sHXHl1H8bjhfP3i6UGXIyKCruAdg1DY+dJDr7O7sY1lnzuHYRl6xF9EgqdAP0ruzneWb+T5t2r5wZWz1NUiInEjqi4XM7vEzLaYWZmZ3dbL/kwzeziyf7WZFce60HjQGQrzrSc38LtV21n8gclcN1cTP4tI/Og30M0sFbgTuBSYCSwys5mHNLsZqHP3KcDPgB/FutAguTtrttdxza9X8vtVFXzuA5P5xqXTMdNdLSISP6LpcpkLlLl7OYCZPQQsBDb1aLMQ+E5k+VHgf8zM3N1jWCsAy16t5NcvbMUBHJzuwO1+B8e733v85V73c6BNz/Ue7SLLOITcaekIMW5EBv997WwWztb95iISf6IJ9AKgssd6FTCvrzbu3mVmDcA4YE/PRma2GFgMUFR0bN0VY0ZkMP2EkWBg3d8ZeX/vOpGTZ8N67OuxHmnU677I99Bj+yn52Xz09AnkZGkURRGJT9EEem/9CoeeeUfTBndfAiwBKCkpOaaz94tm5nPRzPxj+aiISEKL5qJoFTCxx3ohsKOvNmaWBowC9sWiQBERiU40gf4qMNXMJplZBnAtsPyQNsuBGyPLVwF/H4j+cxER6Vu/XS6RPvFbgWeAVOAed99oZt8FSt19OXA38DszK6P7zPzagSxaREQOF9WDRe6+AlhxyLY7eiy3AVfHtjQRETkaGstFRCRBKNBFRBKEAl1EJEEo0EVEEoQFdXehmdUC2wP540cvl0OeehUdk17omBxOx6R3x3NcTnL3vN52BBboQ4mZlbp7SdB1xBMdk8PpmBxOx6R3A3Vc1OUiIpIgFOgiIglCgR6dJUEXEId0TA6nY3I4HZPeDchxUR+6iEiC0Bm6iEiCUKCLiCQIBfoR9Dc5drIys21mtt7MXjez0qDrCYKZ3WNmNWa2oce2sWb2FzN7O/I+JsgaB1sfx+Q7ZlYd+a28bmYLgqxxsJnZRDN7zsw2m9lGM/tSZPuA/FYU6H2IcnLsZHaBu89O4nuM7wUuOWTbbcDf3H0q8LfIejK5l8OPCcDPIr+V2ZGRW5NJF/BVd58BnA18IZIjA/JbUaD37eDk2O7eARyYHFsEd3+Bw2flWgjcF1m+D7hiUIsKWB/HJKm5+053fy2y3ARspnsO5gH5rSjQ+9bb5NgFAdUSbxx41szWRCb+lm757r4Tuv8hA+MDride3Gpmb0S6ZJKqG6onMysG5gCrGaDfigK9b1FNfJ2kznP3M+nujvqCmX0g6IIkbv0SOBmYDewEfhpsOcEws2zgMeDL7t44UH9Hgd63aCbHTkruviPyXgM8Tnf3lMBuM5sAEHmvCbiewLn7bncPuXsYuIsk/K2YWTrdYb7U3f8Q2TwgvxUFet+imRw76ZjZCDPLObAMfATYcORPJY2ek6XfCDwZYC1x4UBoRVxJkv1WzMzonnN5s7v/V49dA/Jb0ZOiRxC5xernvDs59g8CLilwZjaZ7rNy6J6T9oFkPC5m9iBwPt3DoO4Gvg08ASwDioAK4Gp3T5qLhH0ck/Pp7m5xYBvwuQN9x8nAzOYDLwLrgXBk8zfp7keP+W9FgS4ikiDU5SIikiAU6CIiCUKBLiKSIBToIiIJQoEuIpIgFOgiIglCgS4ikiD+P7VswBqGb/h7AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.arange(-2, 20, 0.001)\n",
    "y = F_gamma(x)\n",
    "plt.plot(x, y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Моделирование гамма-распределения"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Задача: получить гамма-распределение с параметрами $\\alpha, \\theta$, предпологая, что есть возможность моделирования стандартного равномерного распределения $R(0,1]$. Пусть $\\chi \\sim R(0,1]$.\n",
    "\n",
    "Будем считать, что $\\alpha=n+\\delta$, где $n \\in \\mathbb{N}, \\delta \\in (0,1)$\n",
    "\n",
    "Заметим, что $\\Gamma(1, 1) = Exp(1)$. Тогда $-ln \\space \\chi \\sim \\Gamma(1, 1)$.\n",
    "По свойству гамма-распределения $\\Gamma(n, 1)=\\sum_{i=1}^{n} \\Gamma(1, 1)$, то есть \n",
    "\n",
    "$-\\sum_{i=1}^{n}ln \\space \\chi_i \\sim \\Gamma(1, 1) \\space \\forall n \\in \\mathbb{N}$, где $\\chi_i \\sim R(0,1]$ - независимые случайные величины\n",
    "\n",
    "Получим сначала распределение $\\Gamma(\\delta, 1)$, где $\\delta \\in (0, 1)$. Воспользуемся алгоритмом [6], основанном на выборке с отклонением. Пусть $v_0 = \\frac{e}{e+\\delta}$.\n",
    "1. Положим m=1.\n",
    "2. Сгенерируем $V_{2m-1}$ и $V_{2m}$ - независимые случайные величины, равномерно распределённые на интервале $(0, 1]$.\n",
    "3. Если $V_{2m-1} < v_0$, то положим $\\xi_m=(\\frac{V_{2m-1}}{v_0})^{1/\\delta},  \\eta_m=V_{2m} \\xi_m^{\\delta -1 }$, иначе $\\xi_m=1-ln\\frac{V_{2m-1}-v_0}{1-v_0},  \\eta_m=V_{2m} e^{-\\xi_m}$\n",
    "4. Если $ \\eta_m >  \\xi_m^{\\delta -1 } e^{-\\xi_m}$, то учевличим $m$ на $1$ и вернёмся к шагу 2, иначе $\\xi=\\xi_m$ примем за реализацию $\\Gamma(\\delta, 1)$\n",
    "\n",
    "Итак, теперь имеем $\\xi-\\sum_{i=1}^{n}ln \\space \\chi_i \\sim \\Gamma(\\alpha, 1)$\n",
    "\n",
    "Тогда по свойству гамма-распеделения $\\theta(\\xi-\\sum_{i=1}^{n}ln \\space \\chi_i) \\sim \\Gamma(\\alpha, \\theta)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[1] https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.stats.hypergeom.html <br>\n",
    "[2] https://kpfu.ru/docs/F1021260618/TViMS.pdf <br>\n",
    "[3] Вадзинский Р.Н. Справочник по вероятностным распределениям. - СПб.: Наука, 2001, 295 с. <br>\n",
    "[4] https://books.google.ru/books?id=c210AgAAQBAJ&pg=PA196&lpg=PA196&dq=%D0%BF%D0%BE%D0%BF%D1%83%D0%BB%D1%8F%D1%86%D0%B8%D1%8F+%D0%B3%D0%B0%D0%BC%D0%BC%D0%B0-%D1%80%D0%B0%D1%81%D0%BF%D1%80%D0%B5%D0%B4%D0%B5%D0%BB%D0%B5%D0%BD%D0%B8%D0%B5&source=bl&ots=lCtxofnVqZ&sig=ACfU3U3lFh_rRICVHSUSt2aFNVvS0jrf-Q&hl=ru&sa=X&ved=2ahUKEwj-oeamq6vlAhVjwqYKHVarBAwQ6AEwBnoECAkQAQ#v=onepage&q=%D0%BF%D0%BE%D0%BF%D1%83%D0%BB%D1%8F%D1%86%D0%B8%D1%8F%20%D0%B3%D0%B0%D0%BC%D0%BC%D0%B0-%D1%80%D0%B0%D1%81%D0%BF%D1%80%D0%B5%D0%B4%D0%B5%D0%BB%D0%B5%D0%BD%D0%B8%D0%B5&f=false <br>\n",
    "[5] https://core.ac.uk/download/pdf/6243682.pdf <br>\n",
    "[6] https://cyberleninka.ru/article/v/generirovanie-rekordov-metodom-vyborki-s-otkloneniem <br>"
   ]
  }
 ],
 "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
}