sem10_likelihood/sem_10_lklh_solution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Дайте мне таблетки от жадности! И побольше, побольше!\n",
    "\n",
    "![](https://cs4.pikabu.ru/post_img/2016/08/12/8/1471008417115067244.jpg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from scipy import stats # more style :)\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.style.use('ggplot')  # стиль для графиков\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Задаченька 1. Функция-няша и функция-бяка!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "а) Найдите минимум фукнции $f(x_1, x_2) = (x_1 - 2)^2 + (x_2 - 4)^2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "      fun: 9.095763293783874e-17\n",
       " hess_inv: array([[ 0.9, -0.2],\n",
       "       [-0.2,  0.6]])\n",
       "      jac: array([ 6.37085673e-09, -2.15944951e-09])\n",
       "  message: 'Optimization terminated successfully.'\n",
       "     nfev: 16\n",
       "      nit: 3\n",
       "     njev: 4\n",
       "   status: 0\n",
       "  success: True\n",
       "        x: array([2.        , 3.99999999])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def fun(x):\n",
    "    \"\"\" какая-то приятная функция \"\"\"\n",
    "    z = (x[0] - 2) ** 2 + (x[1] - 4) ** 2\n",
    "    return z\n",
    "\n",
    "x_init = [0, 0]\n",
    "optim_res = minimize(fun, x_init)\n",
    "optim_res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([2.        , 3.99999999]), 9.095763293783874e-17)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "optim_res.x, optim_res.fun"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Теперь добавим в функцию параметр $a$: $f(x_1, x_2) = (x_1 - 2a)^2 + (x_2 - 4)^2$\n",
    "\n",
    "б) Найдите экстремум функции при $a=4$ и при $a=10$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "      fun: 3.131976557517399e-12\n",
       " hess_inv: array([[ 0.59999985, -0.19999989],\n",
       "       [-0.19999989,  0.90000014]])\n",
       "      jac: array([-1.60278948e-06,  3.16307257e-06])\n",
       "  message: 'Optimization terminated successfully.'\n",
       "     nfev: 20\n",
       "      nit: 4\n",
       "     njev: 5\n",
       "   status: 0\n",
       "  success: True\n",
       "        x: array([7.99999919, 4.00000157])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def fun(x, a=1):\n",
    "    \"\"\" какая-то приятная функция \"\"\"\n",
    "    z = (x[0] - 2 * a) ** 2 + (x[1] - 4) ** 2\n",
    "    return z\n",
    "\n",
    "x_init = [0, 0]\n",
    "optim_res = minimize(fun, x_init, args=4)\n",
    "optim_res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "      fun: 2.3881933437907206e-11\n",
       " hess_inv: array([[ 0.51923086, -0.09615402],\n",
       "       [-0.09615402,  0.98076916]])\n",
       "      jac: array([ 2.14204335e-06, -9.52464809e-06])\n",
       "  message: 'Optimization terminated successfully.'\n",
       "     nfev: 20\n",
       "      nit: 3\n",
       "     njev: 5\n",
       "   status: 0\n",
       "  success: True\n",
       "        x: array([20.00000106,  3.99999523])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "optim_res = minimize(fun, x_init, args=10)\n",
    "optim_res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Но иногда встречаются функции-бяки :)\n",
    "\n",
    "$f(x_1, x_2, x_3) = 0.01 (x_1 - 0.5)^2 + |x_1^2 - x_2| + |x_1^2 - x_3|$\n",
    "\n",
    "в) найдите минимум этой функции устно\n",
    "\n",
    "г*) найдите минимум этой функции с помощью `minimize` (автор задачи тоже не знает решения)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       " final_simplex: (array([[0.03077755, 0.00094606, 0.00094727],\n",
       "       [0.03070671, 0.00094229, 0.00094293],\n",
       "       [0.03070712, 0.00094213, 0.0009429 ],\n",
       "       [0.03087131, 0.00095078, 0.00095323]]), array([0.00220291, 0.00220301, 0.00220318, 0.00220326]))\n",
       "           fun: 0.0022029130300699013\n",
       "       message: 'Optimization terminated successfully.'\n",
       "          nfev: 280\n",
       "           nit: 161\n",
       "        status: 0\n",
       "       success: True\n",
       "             x: array([0.03077755, 0.00094606, 0.00094727])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def byaka(x):\n",
    "    \"\"\" известная неприятная функция \"\"\"\n",
    "    z = 0.01 * (x[0] - 0.5) ** 2 + abs(x[0] ** 2 - x[1]) + abs(x[0] ** 2 - x[2])\n",
    "    return z\n",
    "\n",
    "x_init = [0, 0, 0]\n",
    "optim_res = minimize(byaka, x_init, method='powell')\n",
    "optim_res = minimize(byaka, x_init, method='bfgs')\n",
    "optim_res = minimize(byaka, x_init, method='nelder-mead')\n",
    "optim_res\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "д) нарисуйте хорошую функцию :)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-5.  , -5.  ],\n",
       "       [-4.75, -5.  ],\n",
       "       [-4.5 , -5.  ],\n",
       "       ...,\n",
       "       [ 4.25,  4.75],\n",
       "       [ 4.5 ,  4.75],\n",
       "       [ 4.75,  4.75]])"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.arange(-5, 5, 0.25)\n",
    "y = np.arange(-5, 5, 0.25)\n",
    "x_mesh, y_mesh = np.meshgrid(x, y)\n",
    "xy = np.array([x_mesh, y_mesh]).reshape(2, len(x) * len(y)).T\n",
    "xy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([130.    , 126.5625, 123.25  , ...,   5.625 ,   6.8125,   8.125 ])"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xy = np.array([x_mesh, y_mesh]).reshape(2, len(x) * len(y)).T\n",
    "z = np.apply_along_axis(fun, 1, xy)\n",
    "z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_mesh = z.reshape(x_mesh.shape[0], -1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(40, 40)"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z_mesh.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Линии уровня функции fun')"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEJCAYAAACJwawLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydd3gc1dWH39miuuqr3ntxl+Vu3ACbXgKITggtBEgISb4UWiAkJKQQSINQQksIWVro2MEN27g3yZZkWbJk9d6llbbN98eujVVsy9ZKs9q97/PsI+3uzJ0zZ+/85sy5TZJlGYFAIBBMXlRKGyAQCASCsSGEXCAQCCY5QsgFAoFgkiOEXCAQCCY5QsgFAoFgkiOEXCAQCCY5QsgFAoFTkSQpXpKktZIk9UqSJPo3TwBCyF0ASZIqJUmSR3jdpLRtgtMjSdJLkiRtUNoOF+JBIAKYCUQrbItHIITcdXgKe6U/9hIIJivpwA5Zlg/LstygtDGegBBy10ADdMmy3HDsNXSDoRG6JEnTJUmynBgJSpK0QZKkl4bst8yxb5zj/a2SJFmGbPN7xza3Ot4nOd4vPmGbix2fvTrSCUiSpJIk6YgkSQ8O+dxfkqSuE8reIEnSPyRJ+o0kSS2O716SJMn3hH20ju9rJUkySZJUJEnSDSP449irV5KkHZIkLT/h+8ckSSobhR8zJEn6XJKktiFlJo10no599I59jJIkfQL42T+WXnB8tkeSpGmObVMkSbJJkrRwSBlLHZ+njGTXSOcgSdKrkiR9ccL7qQ4fPnTCZ5WSJD08pJxBv/nQOuH47D7HZ4+dwlfD6twIvpGBc4HbTqwvJzm/L06sTw7bfyFJ0rOO36PRUTfVJzuewI4QctfACxg4w33+DPSN9cCSJGUC956qLEmSvIBngO6TbSPLsg14EbhdkiTphK+uA2yA4YTPrgbCgHOAG4HLsD+RHONJ4E7g+8BU4J/APyVJOnfIYe/D/vSyBGgA3jiZfafgFSAUu/hEA5ePYp8/AXEO+1912L8QKAJygSrgHUmSVLIsHwH+5zifE7kDWOv4/oyRJCkbWAs8J8vyr86mjBPKCgN+wSl+XwejqXPRwFbgTcf/95+hOd8F6oF5wPew14FbzrAMj0MIuWsQDPSMdmNJkq4DMrAL51h5BvgH0HyKbX6A/QL+4DRl/QOIxy6Kx7gD+JcsyycKQBtwtyzLxbIsfwQ8DHzbEb37Yb+AH5Fl+W1ZlktlWX7SceyHGEyn4+nlANACdJzGvpGYCbwsy/JeR1ltp9pYkqQgIN9h3y5Zlt8G3gG2y7L8jCzLxcB3sKcXjvnh70C+Y18kSQoGrgJeOAt7kSQpHbuIvy7L8iNnU8YQfoVdfPec4pijqnMOH5oAo+PpsvMMbdkky/JvHGmZ/wBfACvPsAyPQwi5wkiSpAe02CPK0WzvD/wO+CkjR1DflCSp59gL+OwUZV0OzAVOKgaSJMVgF9DvAdZT2SbLciN2wb3Tse8UYD7DL/4dsiyfWNYW7E8lqUCa4/8vh+yzEZgy5LOXHOfYB1wA3Drk+5QTfeHYdihHgCskSQo81bmdQDKgxh59H8PieAEgy3I99htChuOjD4FO4Fh66CbsN+6hN8aXhtj6IMOJAdZhj3Y/HaXNJ0WSpFnYI97vn2Kb09U5Z7JvyPtaIHKcjznpEUKuPDmOvyWj3P5BoAZ4/STfv489yjz2umOkjSRJ8gaeBh6WZflUUejvgE9lWd44Svuexy6M4dgFfacsy0MvzmHmjPDZ0G5r0gifPcTX52kAPpQkKfSE76sZ7IuZIxznNiAB6HCI55qzsPVk20kAsixbgJf5Or1yB/CqLMumU5zPTOy+HEo2dgH/JfDqGdyATsafgT/Lsnz4FNucrs6NBpnhvtOOsN1Qn8gInTotGqUNEHAp9pxg6Si2TQUeAJbKsiwPTkUfp0uW5RMbyOJG2gj4P+zR1ake7xdjzxlnj8K2Y6zDniO+C7gZexQ3lDmSJKlPiMoXYL+Ay7FfuAPAUuDgCfssGfIeoPHYuUqS9Aj2fOw5fB3pmk/0hWO7QQXIsrxTkiQDcD1wJTAdeOsU51eBPeefAxxyfKbhhGtJkqRoIAQ4URxfBB6UJOluYAb29MxQGof8diPdYPcAdzuOdyF2If7mKew9FTcBKdifZk7GaOrcaGjC/jQBHA8kcrD7UzBGhJArhCRJkdgvyAeA78qjmxj+R8C/ZVneOcbDq7AL7IVDUhxDeQh4Spbl6tEW7LjYX8AeMZqAf4+wWRjwV0mSnsUuJE8AL8qy3AsgSdKfgCckSWrG/qh9DfYbyvlDygmSJCkK8MUeWdsY3Q3xOJIkXYr9N5gvy3KJo+HvVOfXIUnSew77qrGnWq4CvCVJ+j6wGvg19pTNFyfsVyVJ0ufAs8AGWZbPyM4TaHfUFbMkSTcDeyRJ+kCW5fdO2EYjSZLPCe+1jnP1lmX5xEb1h4B7ZFk+VfuMs+rcF8DdkiR9iT2AeAh7Ck3gBMQji3Jci70R5zpZlp8b5T5W4GdOOLYEfCTL8qbTbNeMPbVyprziOMabJxGJd7BfzJuxR7+fAj8+4fuHsEewz2CPwm8CbpJlee2Qcv6C/WmmCLvQ3+hobBwVjh47bwDflGV5tKktsPeWqXPY/03gI+yNhdnYI+Zk4OoRbpIvYBevs2rkHIrjXH8G/N0RGBzj54DxhNcL2PP6h4YUUYC9R9CpcFad+xH2RunV2NttvgTGenMQOJDECkECZyNJUg52Ac6TZXn3kO82AGWyLI+Yu5+MSPa++2myLC87zXb3YO/mFzskMh53JEnSYPd70kQeVzAxiNSKwGk48p6x2FMLG4eKuKciSZIOe2+cHwF/mWgRd3Cs7UHghojUisCZXA+UYc97f1thW1yJvwA7gGIGD3yaMGRZtsqynKnEsQXjj0itCAQCwSRHROQCgUAwyVEqRy4eAwQCgeDsGNaZX7HGzrq6unEpt6iiiYeeX8ejty1lVuapZ4PV6/W0tLSMix1KYewzce/1r3LZtblccUPeGe/vjj4ZK8IngxH+GM5E+SQmJmbEz90utVLbbJ8KIiY8QGFLlKGyrBnZJpOaGaG0KQKBYIJwSyHXalTog/2UNkURyksaAUjJFPMMCQSegtsJeV1zF1FhAahVbndqo6L8UBPRccH467yVNkUgEEwQbqd2dS3dxHpoWkWWZcoPNYq0ikDgYbidkLd09BEe7K+0GYrQ2z1Ad2c/cUmnnPdJIBC4GW4l5MYBMwNmK8EBPqff2A1paugCICJ6rFNUCwSCyYRbCXlHTz8AwR6aH25utAu5PtIzU0sCgafiVkLe2W2fEyg4wPc0W7onzQ32rpfhkSIiFwg8CbcS8vYeIwDBOs9MrTQ3dBEY7IuP70graAkEAnfFrYS8q8cekQf6e2ZqpbW5h7BwndJmCASCCcathNxkti/I4uPlmdOs9/eZ8PPQm5hA4Mm4lZCbrTYAtFq3Oq1RMzBgwdvHM29iAoEn47SrPj8/Xw3sAmoNBsMlzir3TDA7InKtWq3E4RVnoN+ClxBygcDjcGboej/2FVAUw2SxolGrUKmGzfLoEQz0m/H2Fg2dAoGn4RQhz8/PjwMuBl5yRnlni9liQ6vxzLQKgGnAInqsCAQeiLOew58BfgycdCRKfn7+XcBdAAaDAb1e76RDf42Pjw+SJI26bI1GMy52KIUsg6+v75jOyd184gyETwYj/DEcpX0yZiHPz8+/BGgyGAy78/Pzl51sO4PB8ALwguOtPB6TsFssJkxm66gneHe3CfK9vNV0dnSP6ZzczSfOQPhkMMIfw3GHhSUWAZfl5+dXAm8BK/Lz8//phHLPGC+NCovVhqcuKO3trWWg36y0GQKBYIIZc0RuMBh+BvwMwBGR/8hgMNw01nLPBq3G3lvFbLHhpfW8nivePhoG+i1KmyEQCCYYt2oZ/FrIrQpbogzePlr6RUQuEHgcTu10bDAYNgAbnFnmmeDl6LEyYLbi74HzZnn7aOg3CiEXCDwNt4rIA/zsw9O7ewcUtkQZQvU62lp6lDZDIBBMMG4l5McWlGjv7lfYEmXQRwXQ0dqL2STy5AKBJ+GWQn5sgQlPIyIyEFmGliYRlQsEnoR7CbljHvJODxXy8Cj7ghItjpWCBAKBZ+BWQu7no0WjVnlsRB4eZR9Ye2ylIIFA4Bm4lZBLkkRYkC8tHX1Km6IIQSF++PhqqTnaprQpAoFgAnErIQeI0QdQ1+KZEakkSaRkRHCktElpUwQCwQTinkLe3O2xw/RTMiOormgVQ/UFAg/C/YQ8PJB+k4W2LqPSpihCWlYkNptMZVmz0qYIBIIJwu2EPDbc3uDnqemVlIwIAMoPifSKQOApuK2Q1zZ5Zhe8gCBfIqIDOVzcoLQpAoFggnA7IQ8L8kPn60V5bbvSpihGzoxYivfXihGeAoGH4HZCLkkS6fGhlNV4bhe8WfOSGOi3UFRQp7QpAoFgAnA7IQdIiw+jqqGTfg+NSLNnxOLjq2XvtkqlTREIBBOAWwp5enwoNlnmSK1nRuVarZqpufHs23EUm80zu2EKBJ6Emwp5GAClVZ4p5ACz5ibS2d4nuiEKBB6AWwp5sM6HqDAdRZWeK2LT5ySgVqvY/mWZ0qYIBIJxxi2FHGBmehSFZY0eu+ybLsCHWfOT+GrdYcxmz/SBQOApuK2Qz8qMpt9koeRoi9KmKMbSlVn0dPezZ2uF0qYIBIJxxG2FfFpqBBq1ij2H6pU2RTFyZsahjwhgw+pipU0RCATjiNsKua+3lqwkPXs9WMhVKoklK7MoKaijoa5TaXMEAsE44bZCDpCbGc3Rhk6PnZ8cYPF5mahUEhs/L1LaFIFAME64tZDPyY4FYOuBaoUtUY6QMH9mL0hm4+oS+noGlDZHIBCMA24t5HERgSRFB7N5X5XSpijKJfmzMPaZ+OKTg0qbIhAIxgG3FnKAc2YmUlrdSmOb564sn5CiZ8acBNZ8UICxz6S0OQKBwMm4vZAvnh4PwOb9nh2VX3bdbHq7B1j78QGlTREIBE7G7YU8IlRHZkKYxwt5SkYEM+Yk8Nl7+0WuXCBwM9xeyAEWz0yksr6DyvoOpU1RlCtvnENfr4mP396rtCkCgcCJeISQL5mZiFajYvU2z553JDFVz+LzMlnzQSG1VZ678IZA4G54hJAH+nuzaHoCG/ZUYhzw7NXlr7l1Hj6+Wv75/GZkWUxxKxC4Ax4h5ACr5qfRb7Lw5d6jSpuiKIFBvlx1y1xKCuvYttGzn1AEAnfBY4Q8MyGMpOhgPt9W5vGR6NKVWSSnh/PWy1vp6xUNnwLBZMdjhFySJC6Yn0ZlfQeHqlqVNkdRVGoVN3/nHLo7jbz3z51KmyMQCMaIxwg5wJJZifj5aPnwyxKlTVGc5PRwVlw8hbUfH6Rof63S5ggEgjHgUULu663lwgXpbDtYQ21zl9LmKM4135xHVGwQLz2znl7Rt1wgmLR4lJADXLI4A61azfsbxBzd3j5a7vzBCrrajbzx3GalzREIBGeJZqwF5OfnxwOvA1GADXjBYDA8O9Zyx4tgnQ/nzU1hzfZymtp6PO9ONoSUjAguuy6X9/+1i5lzE7nkKr3SJgkEgjPEGTpmAX5oMBiygfnAvfn5+TlOKHfcuPycTGyyzFtrxAhHgIuvmUVqViSvP7eJ5kaRchIIJhtjFnKDwVBvMBj2OP7vBoqB2LGWO55EhOpYMjORjzYepL3bqLQ5iqNWq7jrB8uRbTK/e/R9sVizQDDJkJzZpzo/Pz8J+BKYajAYuoZ8dxdwF4DBYJhtMik7nWp1Ywc3P/ImFy/O5v9uWa6oLa7Cpi+K+O0j77Piwml8/5FLkSRJaZNcAo1Gg8ViUdoMl0H4YzgT5RMvLy+AYRem04Q8Pz9fB2wEfmUwGN47zeZyXV2dU447Ft78XzHvri3g6e+vIjEqWGlzXIL/fVDMmy99yVW3zOWSa2YpbY5LoNfraWlpUdoMl0H4YzgT5ZOYmBgYQcid0taXn5+vBd4F/jUKEXcZvnlpHr4+Gl7/dL/SprgM1922mPlL03j39R3s3HJEaXMEAsEoGLOQ5+fnS8DLQLHBYHh67CZNHEE6X65ZMYU9h+rZV9qgtDkugSRJ3Pa9paRmRfLS0+s5UtqktEkCgeA0OCMiXwTcDKzIz8/f53hd5IRyJ4SLFqYTGerPa5/uw2qzKW2OS6D10vC9h1YRGOzLn365mtZmz10mTyCYDDi1sfMMcIkc+bG81paCKn7/r6+48/JcLlqYobRZinJirq/maBtP/vgDgkL8+NlvLiMw2Fdh65RB5IQHI/wxHLfIkU92Fk6LZ0Z6JP/8vIDWzj6lzXEZ4hJD+f6jF9DW3MMffv6JWCJOIHBRhJBjzwt/+8o8rFaZlz8Sg4ROJGNKNPc9tJLaqnaefvwz+o2evTCHQOCKCCF3EB0WwNXn5rC1sJpdxcqnfVyJabnx3P2jczlS2sSff7Uas0n0IRYIXAkh5CdwxZIs4iICeeGDXfQLsRpE3qIUbvveUor21/LnX60RYi4QuBBCyE9Aq1Fz95V5NLf38daaQqXNcTkWn5vJrfctoXBPNc888TnGPmVH5woEAjtCyIcwJSWClfNS+XDzIYoqRB/qoSxdlc3t319GSUEdTz34EZ3tonFYIFAaIeQjcOvFM4kI8edPhu0YB0Tj3lAWn5vJ9x65gPqaDn714w9oqOtU2iSBwKMRQj4Cvt5a7r92Pk3tvbz84R6lzXFJZuQl8JMnL8XYZ+LJH/9XjAAVCBRECPlJyE4K5xvLslm7q4KvCqqVNsclScmI4KHfXo63t5bfPvQRhburlDZJIPBIhJCfguvOn0ZaXCh/e28Hze29SpvjkkTFBvPw768gMiaIZ59Yzea1h5Q2SSDwOISQnwKNWsUPrl+A1Srz+ze/wmwRCy6MRFCIHz/99WVkTo3m5Wc2YHh1G1armLdGIJgohJCfhmh9AN+9Zh6lVa38Q4z6PCm+fl488PMLWXZBNp+9u5/fPfwxHW3iKUYgmAiEkI+ChdPjuXJpFp9vK2PdLjFH98nQaNV8894l3PHAcipKm/n5/e9SUihGyQoE440Q8lFy46rpTEuN4Pn3d1Fe26a0OS7NohUZPPL0lfj5e/Pbhz/m47f3YrMpMsumQOARCCEfJWq1ih/esJAgfx+eemMzXb1iJsBTEZcYyqNPX8mcRSm8+/oOnn3ic3q6+pU2SyBwS4SQnwFBOh9+cstiOrr7efrNr7CIBr1T4uvnxd3/dy433b2Ig/tqeOz771J6sF5pswQCt0MI+RmSFhfK3Vfmsb+skRc/2I1CC3NMGiRJ4tyLp/LgU5ejUqv4zYMf8c7r27GYRQ8ggcBZCCE/C1bkpfCNZdms2V7Oh5tEv+nRkJIRwePPXsXiczP45O19/OIH71FxuFlpswQCt0AI+Vly46rpLJgWz2uf7hMjP0eJr58Xt31vGfc/sorurn6e+NH7GF7ZhmlATIkrEIwFjdIGTFZUKon7r51He5eRZ/6zlSCdN1NSIpQ2a1Iwc24SGVOiMbyyjc/e28/urRV867tLyZoWo7RpAsGkRETkY8Bbq+HBW88hIkTHk69t4mhDh9ImTRr8/L259b6l/N8vL0GW4akHP+L1v20Sc5wLBGeBEPIxEuDnzaO3L8XHS8PjL22gvqVbaZMmFTkzYnniz1ez6orpbFhdzIPf+Q/bNpaJRmSB4AwQQu4EIkL8efT2pVisNn7+4nqaxARbZ4S3j5brbl/Aw7+7guBQf/7++7X89qGPOVreorRpAsGkQAi5k0iMCuaxO5bTN2Dm0RfW0dopVs45U1IyInjk91dwyz3nUHO0jccfeJcXn15Ha3OP0qYJBC6Nywj5vromylraaeqZvAKYEhvCo7cto6t3gJ+/uJ6ObjGS8UxRqVUsvzCHp164jguvmsmOzUf42d1v8c7r20X+XCA4CZJCuUi5rs4+mVJ1Rzd/3boXrVpFqK8PpS3tPHruAqIDdeNuhF6vp6XF+Y/vRRXN/OLlDUSG6XjirhUE+ns7/RjjxXj55GxpaermvTd2snXDYQKCfLj8+tksXZWNRqOeMBtczSdKI/wxnInySUxMDIA09HNFI/KSplZ+/+UObLLMEysX88A5eSxLiecPm3Zhtk7ekX85yeE8eOsSGlp6ePylDXT3iXlZzhZ9RAB3/XAFP//jN4hNCOWfz2/hkfveZueWI2IiLoHAgWIR+YGyI/z40w2E+/uRGR5KTWc3P1s+D7VKxaHmNjLDQ8fdiPG+i+45VM+vX9tEbHgAj96+jNBA33E7lrNw5WhLlmX27zzK269up666g7ikUC66aiZzz0lFrR6/mMSVfaIEwh/DGU+fdLb3sebDQoy9Jn76y2vBlSLy9w8eZnFSHL++cAm35k2lx2SioME+ZPtQcxuflhzhhe37lTLPKeRmRvPwt5bQ2NbLg899QX2r6Jo4FiRJYubcJJ748zXc+YPl2KwyL/xhHT+96y3WfXJQjBAVuBVNDV28/rdN/Oj2N/ns3X30nWLGVcUi8kfe+5TZsZHMjY+mobuXdwoPcUl2Kq/tPojRbGFRYgx765qwyTKPnb9oXIyYqMiitLqVX/5jI2qVikdvX0pyTMi4H/NsmUzRls1mj9A/eXsv5YeaCAz25fzLprHiwhz8dM5rl5hMPpkIhD+G40yfVFe08um7+9i+qRyVSmLRigwuvGomUTFBJ82RKybk723bxcs7C1meEo9WrWZvXROljpTKry44B4Cu/gGeXL+Nny2fT5CP8xsMJ7JCVjd28vjLGzD2W3jw1nNcdjj/ZLxIZVnm0IF6Pn1nH4V7qvH182L5hTmsvHwaQSF+Yy5/MvpkPBH+GI4zfFJ60F6H9++qwsdXy9JV2ay6YjohYf7HtzmZkCs218r8BPu8GmUtHdhkmbY+I+emJ/Kd+TOPb1Pa0k6kzh8/7eSfEiY+Mojf3HMej720kV+8vJEf3riQuTmxSpvlFkiSRNa0GLKmxXC0vIVP393HZ+/vZ80HBcxbksZ5l04lKS1caTMFgmFYzFZ2bjnCFx8f4MihJnSBPlx50xxWXJSDLsBn1OUo3v0QwGg288quA1w9LZMInT2CKmtp55E1m7kldwoXZqWMixFKRBZdvQP88pWNlNe2c9fls1k1P21Cj3863CXaaqzrZM0HBWxZV8pAv4Xk9HCWrspm3pI0fHy1Z1SWu/jEWQh/DOdMfdJY18nGNcVsWVtKV4eRyJggzrtkKuecn4m3z8nrp8ulVo4JuSzLlLd28JsN2/nz5efhq9Wwu6aBJ9Zu5VtzpnF5Ttrx7SRpmP1jQqkKaRww84c3v2J3ST2XL8nilgtnoFI599zOFne7SPt6B9iytpQNnxdTV92Oj6+W+UvTWLoqe9RRurv5ZKwIfwxnND4xm63s2VbBxs+LKS6oQ6WSmDk3kWUXZDNlVvyoNMBlhfwYz23dy/76ZnIiwyhuauPaGVksS4kHwGqzoVZ93cHGWaKuZIW0Wm28/NFePtt6mPlT4/j+tfPx9lI+heSuF6ksy5SVNLLx82J2bC7HbLKSmKpnycos5i5ORRd48sdYd/XJ2SL8MZxT+aS6opUt60rZsq6Unq5+9BEBLFmZxeLzMgflv0fDuAp5fn7+BcCzgBp4yWAw/OY0uwwTcoBtVXVE6vzx1qiJCdRhk2VkWUatUmG2WukwDrC27Chmm42siDDmxEWNyW6lK6Qsy3y8pZRXPt5LcnQIP71lMeEhZ/bDOhulfTIR9PYMsHXDYTauLqamsg21WsW02fHMX5bGzDmJwx5tPcEnZ4Lwx3CG+qS1qZttG8vYurGM2qP2OnYs+s6ZGXfWT+DjJuT5+flqoBQ4H6gBdgLXGwyGolPsNqKQH2NoBN5rMrPxSDUVbZ3UdfWwPDWeP3+1l99etHRMA4dcpULuKq7jj//eilar4sc3LSYnWbmGOVfxyUQgyzJVR1rZuuEw278so6OtDx9fLbkLklm4LJ3s6TGo1CqP8sloEP4Yjl6vp7Kihl1bjrB1w2FKDzYAkJYdyYKl6eQtTiEwaOwDAsez18pcoMxgMBwByM/Pfwu4HDiVkJ+SgoZmqjq6uTwnjX6LhfXlVVR1dJGuD+bOedPxUqs52t5FTWf3hIwAHW/ysmN46r7zefK1Tfz8xfXceflsVs5LVdost0eSJBJT9SSm6sm/dR4lB+rZtvEwu7ZU8NW6UoJC/Jh3TirLL5hJeIzPqEaPFhfUUl7SSFxSGGlZkegCfcalfUfgGvR093NgTw37dqxn11dlWC02ouOC+cZNc5i/NI3wqMAJscMZQh4LnLhoZQ0wb+hG+fn5dwF3ARgMBvR6/UkLPDcsjP01Dej1epp7eqns7mNuSiIXTckAoLCukW01jayann3Kck6HRqMZ0/7ORK/X8/KjMTz2whqee28nRxt7eODGJfh4n1kPi7HiSj6ZaCIiI1hy7gxMAxZ2flXGxtUHWPdZEWs+LEQX4EPu/BS+dd+56CNGvjgtFitqVQvGXhtb15ezdX05t967gvikcI6WN1FT1UZKeiTRcSGTWtw9uY7IskxVRQu7vipj55bDFBfWYLPKhOp1XHrNHJaumkJqRtSE/7bOEPKRLB6WrzEYDC8ALxz7/nSPZvG+WlpaWnh+2z76B0zMjQylpaWF0pY2NlXUcm5KHHE+Wqrr67HJoFWp8DrDGfFc8RHxxzfO5z9fBPD22oMUHannxzctIlofMGHHd0WfKEHmND2Z05Zxc99Cqsq72Ly2kMK9RzH299DSMvJ0urIsk5IVQsZUu8g9+8vP2buzlPq6JjZ8XkxHay+tzd1cdcs8cucnTVox97Q6YjZZKCmsY9/OKgp2VtHSZJ9qIyEljIuvnsWMvATyFmTT3t4GQGtr67jZ4kitDMMZQl4DxJ/wPg44eQL8DLDJMhIcT59sq6qjsKEFrUrFjJgIihpb+OPm3aSGBnO0o5NnLz33jMXc1VCrVNywchqZCWE8859t/OhPa/hu/jzmT41T2jSPxNfPi0XLs8icpj+t8EqShEajxmy2otWq0WrV+Af48O7rO1i0IoOFKzI4uLeGvTsqyZwajb9jGt4goZYAACAASURBVIFj5TY1dLFl7SFsNpmlK7PRR07cDVzwNTarjcryFooL6igprKX0YAOmAQteXhpyZsZy8TUzmZ6XQKj+66m2x3PSttHgDCHfCaTn5+cnA7XAdcANTigXlSSxKiOJn322id21jTT29LIyPZmsiFAae/r478HDXJadyqU5aTy3bR8v7ijg3oWznHFoxZmdFcMfvreK3/1zC0+9sZnLzsnkpgumo53kN6rJzGijZ61WzfYv7b0VktPCqapo5YfL0rFZbaRkRvDmi19x7be+/h2PlfvCH9YxPS+BnrZePjLs4brbF+Dr54XNamPbxjK03hqyp8ccH/E3WSN6V8Nmk6k52kZJQS3FBXUcOliPsdf+1BUTH8Li8zKZkZdA1rQYvLyV7yI8EmO2ymAwWPLz8+8DVmPvfvgPg8FwcMyWOUgODeaZy1bQ3NtHckgQOm8vNhyp5uPicm6clXO8C6LOS4vJanPWYV2CiBB/nvzOubzy8T4+3HSIoopmfnD9gglNtQjOjJ7ufla/X8Dh4gZ+9MQltLX04Ofv5ehuJtHZ3ke/0YzWS4PNJqNSSdisNr7acBh9hI7Lrs0F4KfffgvTgAVfPy/eeH4zJpOVnq5+dmwq584HluPlrUGSJMoPNeLn50V0vOtOxOZqDPSbqTjczJHSJsoPNVJ6oIEex2peEdGBzF2cQvb0WLKmxThlrp6JwCm3F4PB8CnwqTPKGomoAH+iAuz9q/fXN/H8tn38dNk8ZsZEIMsyNZ09tBv7WZocf5qSJh9ajZq7rpjNtLQI/vrODn7w7GruuCyXFXnJIhpzMSrLmtmythT/AG/u+9lKdIE+7NxcTnL6191Jy4obiYyxN5babDZUKjXdXf001HaQPd0+905HWy9Tc+OorWqjq9NISWEdv37+OgCeefwzqitbSc2MZM+2Sv73USH9fSb8/L25/5ELXDZiVAqL2Up1ZRuVZc1UHG6m8nATtVXtxxclCY8KZMacBLJn2IU7LHz8VyYbDybdr240W7gkK5WZMREMWCw09Rj5oOgwsgwzYiLYU9uIyWql3djPhZnjM0eLEiyYGk9abCjPGrbxl3d2sL2olnu+MYfgM5hYRzB+9PUM8Pffr6O3p5+rb5mHsc+ELtAHbx8tQcFfR3VlJQ3HBVt2iElnhxGzyUqoQ0S6u/rRatUM9FtobW4mPScasA9kypwWzdHyFoJC/Fj78QGuu20Bial6DK9sY/fWChYsS5/gM3cNbFYbTQ1d1Fa1U1fV7vjbRn1NBxaL/UldF+BDUrqemfOSSM2MIDkjwil9u12BSSfkJquNzw4dYV5CNAcbWzja3oVWreK+xbk8t20fpc1tTI8OZ0d1Aw3dvXwrb5rSJjuN8BB/fnHnCj7afIh/rS7g/j9+xne+MUc0hLoAPn5e3Hb/UmqPtlFW0sj6z4pYdkE2s+YlcXBfDZ++uw+zyYpKpSJvYTIAWseUDFqNCrPJiq+fvatpa1M3kiShC/DhyKEmklLtvWBUKonuzn5C9ToOFzXg46cl0fFdeFQgRftrWbAs3W1z57Is09VhpLmxi+aGbpobu6mraqeuut0u2Oavl4cMi9ARmxDK1NnxJKeFk5QWjj4ywC39ApNQyJckx9HY3curuw4QHxzAjJgIlqXE86t1W6nr6uHpS5bjrdFwXloSf9y0i67+AQK8vdzmB1SpJC5fkkVuZjTP/GcbT72xmWW5SdxxWS7+vl5Km+exqFQS6dlRpGdHsewC+2fHBHXRuZns/qoCs8nCtbfNJzjUn9bmHspKGsidl0R0fAiHixs4/zJ70LFn21Gi44JJSAnjk3f2kjXd3uXM2GeivbWXRSsyWPNhISnpX89pX13RSkDgsUZQmGzVXZZl+o1mOjuMdLX30dneR1tLL82NXbQ4RLulsRuTafAqUGHhOmISQpgyM5bYhFBiEkKIjgvG18+zroVJJ+QA10zP5OKsFPy87BHMs5t30z1g5q9XnH98m6aeXtLDQwgchwUpXIH4yCCeuvd83l57kHfWF1FY3sh9V89lZka00qYJHBwLHmbkJTAjL2HQdxqNigGj5XhUfs75Wbzx3CZiE0Opq27nyhvz8PbR0trcQ7Cjwe1wUQNqtYrI2CDqqttZuPzrNMrRIy1cfNVMx3EH2/H8776gq8NIUIgfukAf/Py88PHzwtdPi6+/F76+Xvj6eeHjp8XXzwutVo1arUKtUaFWq1CpJdRq1YjBkM0mY7PasFptWC2Ov1YbA/0WjH0m+o1mjH0mjL0mjEYTxr6v33d29NHZfky4jcNEGuzdP/WRAUTFBjE1N47wqEDCIwPQRwaij9CdcspXT2JSCjlwXMTb+oz0Wyz8eOnc49819fTx1617j0+B665o1CquXzmNvOwYnv3Pdh5/eSMrZidz6yUzCfBzzxuYuxAU4seSlVnHo/blF+YQERVAQ20nKy7MISTMH1mWWXXFdN7/1y6mzY5n25dlXHzVTDQaNcGh/lgdvbR6ewbo6xkgNSsSGN5N0s/fm9bmHo6UNtHd2U+/0cTZTLF0oqhbHKJ9LM9/RuVoVPj6eREU7EtQiB+p2ZEEBfsRFOLr+OtHYIgvIWH++Ou83eZpejyZtEJ+jIbuXuq6egn1sz9Wlrd28NPPvuSS7BQun+IZDT/p8WE8ff8qDGsP8v7GYnaV1PHNi2eyPDdJXAQuzrHfR6tVM3Nu0rDvZi9Ipruzn4rDTSw5P4upufaeWQuXp7P+s2I0GjX7dh5lel7CSbvK3XLPOYPe22wyA/3mr6PlPhP9jki5r8+ExWzFZpWxWm3HBdt2QsTt7e2DyTRgF3WNCs2g6N3+19tbY4/2/bzw8bVH+r6OpwCtC0zX7G64zHzkY+HBz79EJUkkhwaz9Wgtl2SncoVDxE/V8OOOQ40r6tr5+/u7OFTVSk5yOHdfmUd8ZNCo93dHn4wVV/RJT1c/H7+9l4rDTaRlRXHlTXloJmiwmCv6Q2kmyicuv7DEmAqTZT4oKiM2UEeAtxdZEWHHPz9VROquFdJmk/li5xHe+Gw/xgEzly/JIv/cKaNauMJdfTIWhE8GI/wxHKWF3C2ecSRJOh6BD/3cE1GpJFbOS2XelFhe+3Qf720oZtP+o9x52WzmiAWfBQK3Q9mZXiaQirZOzFbr6Td0I4J0Pnwvfz6//PYKfLQannxtE798ZSO1zV1KmyYQCJyIRwh594CJBz5ax08+/ZLO/gGlzZlwpqRE8If7V3HrxTMprmjh/qc/4x8f7aW7z/N8IRC4Ix4h5AHeXty7cBYHG1u45/3/caStQ2mTJhytRs3lS7L46/9dzPLZyXyypZTv/PZj/ruxBJPZs55UBAJ3wyOEHOD89CSeuWwFJquV+/67lk0VNUqbpAjBAT7ce/Vcnr5/FZkJel77dB/3/eETNu6tPD6RkEAgmFy4Ra+VM6Gl18ijazZT0tzGd86Zw1VZyag8tFEUoKCsgdc+2c+RunZSYkP47nVLSAz38diG4pEQvTQGI/wxHKV7rXickAOYLFae3rSLNYcrmZ8Qw4PL56Hz9qy5GU7EZpPZtP8o//q8gOaOPqamRHDDqmlkJ4WffmcPQAjXYIQ/hiOEXCFkWeZ/lfX87ovNRAb48/j5i0gNC1bUJqUxma18VdTAax/upKOnn1kZ0dywahppcaFKm6YoQrgGI/wxHKWFXP3YY4+N+8FH4LHu7m4ljnscSZKYl55MVoiOtWVH+e/BMkL9fEgLC/bYtIJarSJvagpLZsTg76NlS2EVn2wppby2jagwHWFBk2O1FGfj5+dHX1+f0ma4DMIfw5konwQEBAA8PvRzjxVysDtfp5I4Ly2RkuY23jtwmLquHmbHRqJVe+bamH5+fpgG+slOCueC+Wn4eGnYUljNp1sOU1rVQmSoDn2wZwm6EK7BCH8MR2kh99jUCgx+HLLabPxrbzGv7zlITKCOn5+3gNQwz1sHcaRHROOAmc+2lvHBlyV09Q6QkxzOFUuzmJ0Z41iL0r0RqYTBCH8MR+nUihDyIc7fV9fEr9ZtpWvAxL0LZnFpdqpHpVpOVSH7TRb+t72cjzYformjj7iIQK5YksWSWYloJ2jCJiUQwjUY4Y/hCCFXkJM5v93Yz2/Wb2dnTQNLkuP44ZI5BHhIr5bRVEiL1caWgir+u7GEyvoOQgJ8uHRxJivnpbrlKkVCuAYj/DEcIeQKcirn22SZ/+wv4R87Cwn18+Eny+aRGxs5wRZOPGdSIWVZZv/hBt7fWEJBWSO+3hpWzkvjooXpRIT4j7OlE4cQrsEIfwxHCLmCjMb5pc1t/HLdNmo6u7lySjp3zpuOj8YtJo0ckbOtkOW1bfx3YwlfFVaDDLOzolk1P42ZGVGoVZN7ALEQrsEIfwxHCLmCjNb5/RYLL+0o4L0Dh4kL0vGTZfOYEqmfAAsnnrFWyOb2XtbsKOeLHUfo6OknIsSflfNSOTcvheAAHydaOnEI4RqM8MdwhJAryJk6f29dI7/dsIPmXiP50zO5NW8qXm7WTdFZFdJssbKjqJbPt5VxoLwJjVrFvCmxXDA/jSkpEZOqAVkI12CEP4YjhFxBzsb5vSYzz2/bxyclR0gOCeIny+eSoXefkY/jUSFrmrpYs72c9bsr6DGaiA0P5Lw5KSzNTSQkwNepxxoPhHANRvhjOELIFWQszt9WVccfvtxJh3GA62dmc1NujltE5+NZIQfMFrYUVLN6WxmlVa2oVBKzM6NZnpdMXlaMy3ZhFMI1GOGP4QghV5CxOr+rf4C/bd3HmsOVJAQH8MDiPGbERDjRwolnoipkdWMn63dXsGFPJe3d/QT4ebFoRgLLZiWRkRDmUqkXIVyDEf4YjhByBXGW83dU1/PM5t00dPeyKiOJu+fPJMjH2wkWTjwTfZFarTb2ljawcU8lO4pqMVmsRIfpWDIrkUXTE4iPDJowW06GEK7BCH8Mx2OF/BebPuW2lBmKzgXuTOf3Wyy8sacIw/4S/L20fHv+DC7ISHapyHI0KHmR9hpNbDtQw8a9lRw40oQsQ3xEIAunx7NwWgIJUcqIuhCuwQh/DMdjhTzp9d+wSB/Hz6cuRqdRZjTgeDi/oq2DP27azYHGFqZHh/PA4jwSQwKdeozxxFUu0rYuI9sOVLOloJriymZkGeIiAlk4LZ6F0+NJiAyasJukq/jEVRD+GI7HCvmftq/j2dKdxPjq+PX05aToJn4u8PFyvk2W+fxQBX/fvh+j2UL+9ExunJWDr9b1BxK54kVqF/UaviqsoqjCLuqx4YHMnxrH3JxY0uJCx3XyLlf0iZIIfwxnPH1iHDCz91ADJouVGy5ZCK4k5HV1dexrb+Thwo0YrRYeylnIisikCTVivCtkh7Gf57ftZ83hSsL9fbl7/kyWpcS7dLrF1S/S9m4j2w/UsKWwmqKKZmw2mWCdD3nZMeRlxzA9LRJfb61Tj+nqPplohD+G42yfNLX1sKuknl3FtRSWN2Gx2kiJCeHfT30LXE3IAZr7+3iwcAMHO1vIj8/m3vRctKqJ6YY2URXyQEMzf9qyh7LWDqZEhnHPgllkR4SN+3HPhsl0kXb3DbDnUD27iurYc6ievgEzGrWKKcnh5GbFMDsrmhh9wJhvnJPJJxOB8MdwxuoTs8VKydEWdhfXsftQPTVNXQDE6AOYkxPLnOwYshL1xMfHgSsKOYDZZuXPh3fzTnUJUwL1PDF9CVE+unE3YiIrpNVmY3VpJS/vLKTd2M95aYncMXc6ETrXWqRhsl6kZouVQ0db2FVSz56SOqodF4I+yI/paZFMS4tkelokoYFnPgBpsvpkvBD+GM6Z+sRmk6ms76CgrJGCsgaKKpoZMFvtgUhKOLMzY8jNiiY2fHD7msvlyEfqfriusZIni7aiUUk8nLOIxeHx42qEEhWyz2Tm3/uKMRQeQoVE/oxMrpuRha/WuemAs8VdLtKmth72lDZQUNZIYXkjPX0mwN4L5piwT02JGNW0u+7iE2ch/DGc0/lElmUaWnsoKGtkf1kjB8ob6XbUyThHnTz2OlVqcFyEPD8//3fApYAJKAe+ZTAYOkax60n7kVf3dfFwwUYO97STH5/NPem5eI1TqkXJCtnQ3cuLOwpYX15FmJ8v38qbyqqMJMVnCnTHi9Rmk6mob7eLelnj8ehHJUmkxIaQnRROTnI4WUl6gnXDJ/ZyR5+MBeGP4Qz1ic0mU9PURXFlM0WVzRRXNNPcYV8KLizI1yHaUWf8lDheQr4SWGcwGCz5+flPARgMhp+MYtdTDggasFr5W9lu3q4uISMglF9MXUKCv/O78LlChTzY2MJzW/dR1NRKYnAgd8ydzsLEGMUaRF3BJ+ON2WLlUFUrBWWNFFU0cbiqDZPFCthzktlJerKTwslODic6TEd4eLjb++RM8IQ6cqYEBoWwfX8pxZUtFFc2U1LZQo/RHnGHBPiQlRTO9NQIpqVFjqndZtxTK/n5+VcCVxsMhhtHsfmoRnZuaq7mV0VbMNtsPJA5l4ujnbvsmqtUSFmW2VxZy0s7Cqju7CYnMozb50xjVszEL2ThKj6ZSMwWK+W17RRXNh+/EI+lYoJ1PuSkRhGn9yc1LpTU2BBCA33HVA9tNpl//6+QjPgw5uTEYrPJtHb2EeDvjY+X6KLq6litNqqbujhS20Z5bTvlNW1U1HUcDwZiwwPtwUByONlJeqJCdU7TrYkQ8o+A/xgMhn+e5Pu7gLsADAbDbJPJNKpyG/q6+b9tn7G9qYaVcWk8Med8QrydM2OeRqPBYrE4pSxnYLHZ+O/+Yp7fvJOm7l7mJcVx39J5zIiNmjAbXM0nSmCzyVQ1tFNwuJ7CsjoOVTZztL4dm+NaCQ30IyMxnMzEcDISI8hMDCdilBfrgMnCF9tLefmD7cydmshPb13B6q2HeOOTXQAsnpnM3VcvRJblYeW1dxuxWm2EBvopuui1J9URs8XK0fp2DlU2cehoM4eONlFW3YLJbBdtX28tGQnh5KRGMTU1imnp0eM6o6eXlxecjZDn5+d/AYykJA8ZDIYPHNs8BOQB3zAYDKO5M5zRXCs2WeatqiKeL9tLsNabh6YsYl5YzKj3PxmuGlkMWCx8VFzOm3uL6egfYH5CNN/Km0a6PmTcj+2qPlESvV5PTV0DFXXtHKltp9wRidU0dh0Xdz9vLQlRQdx5+WxSYk/+O+0rred/O46QFBOMRqUiLzuGFz7YzY9vWoTVKvPqJ3u5aGEGGQmDu6fWNnfxz88LqKhrJzzEn3uvmktUmL1nV0VdO5v3V5EUHUxuZvRJG3BHujmcrT/crY5YbTYa23qpauikqqGDqsZOqho6qWvpxmqz/8a+3hqSY0JIjQ0lNc7+N0YfgEolTf6Rnfn5+d8E7gbONRgMfaPc7awmzTrc3cZjBzZT0dvB1fFZ3JOWi4/67B9FXb1CGs1m3j9wmLf2l9BjMrMgIYabc3PIGsc+6K7uEyU4mU8GTBYq6juorGvnqOPCv++auUSHBYxYzr7SenaX1JOXHcOR2nZ8vDX4emvYfrCWn9y8mAGThQ83HQIJrlkxBavNhlqlwjhg5r8bS+g1mrjj8tl8sqWUksoWfnjjQoorm3lvQzGx4YE0tPaQlajniqVZ2GwyBWUNvPbpfoL8vblgQTrzp8aNqz8mA8YBM3Ut3dS39FDX3EV9aw/VjZ1UN3Udj7IBIkP9SYgKJiEyiMSoIFJiQ4gOCzjpk5DSQj6mhFx+fv4FwE+ApWcg4mdNekAo/5h7Ec+V7cVQXcy2lloembKIacGTe+rYk+Gr1XLDrBwuy0nj/YOHeaewlHv++wVz46O4OXeK2y43N1nw9tKQlagnK/H0v0Ov0URRRTMLpsWRkxzB3tIGEqODqWnqIi4i8Hh5fQNmvLX2XlrHYqwjte0YB8wsmpEAQICfFz7eGswWKwVljeiD/bj14pkUlDXy8ZZSrliaRWF5I+t2V3DHZbk0tveyu6SOnORwAv0n56yco8Vmk+no6ae5vZeWjj4a23upb+l2iHc37d39g7YPC/IlPiKIC+ankRAZREJUEHERgU4fHTzejLVl5S+AN/C//Px8gG0Gg+HuMVt1CrzVGr6fOYdzwuP5VdEW7t71Odcn5nBnyky8xxCduzI6by9uzp3CN6Zm8MHBwxgKDvHdD9aSGxvJzbNymB4d7tLD/gXQ0NbDtgM1bDtYS1SYjqqGDny9NRw62sLCaV+Pl+jo7ic7KXzQvm1dRtRq1fHcq02WCdL5UF7TTr/JwpRkeyAT4OdFRIg/tc1dHDragq+3likpEST2mSirbqOgrJHFMxKclmKZaMwWK509/XT0DNDR3U9Ht5Hmjj5aOvpo6rALd0tHHxarbdB+wTofovU6cjOjidYHEON4RYXp8J4EjcujYUxnYTAY0pxlyJkyOzSKN+Zfxl8O7+LNo0Vsbq7hZzkLmBE88T09Jgp/L3uEfuXUdD4sKsdQUMIDH69nSmQY18/MZkGCct0WBacmJSaEx+9aTktHH4er22ho7UGtkmjp6EPnZ4+SWzv7GDBbidHbc9/H0p79JgteGjVajX2MQXefCV9vDSaLBbPFRohjUWtZBovVhslspb27n0THtL86Py/MFit9/eaJPu1h2GwyZouVfpOF3n4zvUaT42Wmt99Ej9HxWb+Zrt4Bu3B399PR0z+i/ZIEoQG+6IP9SIsLZcHUOMJD/NEH+xER7E94iD9+PpMruj4bJvXtyF+j5SfZC1gekchvirfynV2ruSouk7vTcvHXuO+P56vVcu2MLC6fksZnJRX8p6CEh1dvJjEkkKunZXB+WhJeLrpsmqciSRIhAb6EBPgSrPNha2E1Vy7NJjJUx8EjTWQl6vlixxGiw3QkRNlnAj229F20PoDiyubjN+mqhk4yE8II1vnSazQdF6q+fhOyTSbAz5vWzj5ys6KPH7+9u5+QwOGDncCe9vnp377AS6vGS6tGo1ahliRUqmMvFSrHe7Uk4ePrjdHYj80mY5NlrDbZ/v/x9zbMFvsNZcBsxWSyMGCxYjJZj3fROxVqlYS/rxcBfl4E6XxIig4mOMCHYJ0PQTofgnTeBOt8CAnwISTQ12WXCJxIJrWQH2NuWAz/nH8Zfy/fyzvVJWxuqeEn2fOZHxartGnjio9Gw5VT07k0J5X15VW8XXCIP3y5i5d3FHL5lDQuy0kjxHfki1egDLIs4+2l4aJF6VhsNhbPSGBvaT2PvLCOqDAd91w1h0B/b0qrWqlu7OScmYnkJIXz1v8KqWroQJKguLKZSxdnkhAVxOHq1uPpgb2lDYQG2qPTpvbe46NUu3oHsFisRAT7A4z41BYfGeQQXos9arZascnyiGKtUqmQZdvX4q5S2QX/BPH39dYQrPPB20uNl0aNt5cGb60ab60GL60aby81fj5e+Ptq0fl64X/C/15atXiyPENcaq4VZ1DY0cSTRVs52tfJyqhkvpeRR6jXyP06J3Pr+0jIssy++ibeLihlW1UdWrWK89OTuHpaBkkho1tdx9184gyU8El7t5H6lh5yku358p3Ftby99iB9/RZuvGAa86fEIUkS764voqqxk6ToYLYdqOG+q+cSHxnEX9/ZQXZyOCtmJ/PpV6VUNXbxzYtmOKURT9SR4Sjda8XthBzsQ/xfryzkjcoD+Ko13Js+m0ti0oYtK+fOFbKqo4t3C0tZXVqJyWplbnwUV03LJC828pTRjjv75GxxZZ8YB8x8seMIzR19zM6KZka6fchHc3svT72xGRnQ+Xpx1xWzh82kd7a4sj+UQgj5OFLZ28Fvi7ezr6ORGcER/ChrHqm6rwdreEKF7Owf4KOicv578DBtxn7igwK4NCeVlelJBI6wQLQn+ORMmaw+6TdZaO3sQ5IkYvQj920/GyarP8YTIeTjfSBZ5pO6Mv5atocei4mr47K4PXUGOo2XR1VIk9XKhvJqPig6THFTG15qNctT47kkO5WciLDjUbon+WS0CJ8MRvhjOEoLuVs0dp4KSZK4JDadc8LjeaF8H4bqYr5orOTe9NncGOaaq/SMB15qNSszkliZkURZSzsfFZfzRdlRVpdWkhIaxKU5aZyXlogYYiQQTD7cPiIfSnFXC38o2UFRVwt54bHclzKLjIBQRWxRmj6TmXXlVXxYVEZZawc+Gg0XTklneWI0UyL1oueAAxGBDkb4YzhKR+QeJ+RgHxn3cV0Zfz+yj44BI5fGpnNX6syT9m5xd2RZpqS5jY+Ly9lwpBqj2UJcUAAXZCRxfkYS4f6utRzdRCOEazDCH8MRQq4gXoE6frdrPe9Ul+Cj0vCtlOlcE581YYs/uyJ+gYG8u3Mfnx+qpLChGZUkkRcXxQWZySxMjMFL7Xm+EcI1GOGP4QghV5Bjzj/a28mfSnextbWWWN8A7knLZVlEgkemFk6skLWd3XxeWsnq0gpaeo0EeHuxLCWec9MSmRqlH9ad010RwjUY4Y/hCCFXkKHO39ZSy58P76ait4NpQeF8NyOPqUHhpyjB/RipQlptNvbUNrLmcCWbK2oZsFrR+/uyLCWe5akJZIWHuvVNTwjXYIQ/hiOEXEFGcr7FZuPT+nJeLN9Hq8nI8ohE7kqdSaL/6EZGTnZOVyGNZjNfHa1jfXkVO6obsNhsRAf4syw1geUp8aSGBbudqAvhGozwx3CEkCvIqZzfZzHz5tGDvFVVxIDNykXRqdyeMoMIH/8JtnJiOZMK2TNgYnNlLevLq9hd24hNlokPCmBZajznJMW5jagL4RqM8MdwhJAryGic32Yy8lpFIe/XlKKS4Or4LG5OnEqQl3tORnW2FbLD2M+mihrWl1ezv74JGYgO8GdxchznJMWSEzl5c+pCuAYj/DEcIeQKcibOrzf28NKRfXxefwQ/tZYbEqeQn5DtdtPlOqNCthv7+epoLZsqatlT24jFZiPU14dFSbEsSoplZkzEpOr9IoRrMMIfwxFCriBn4/zyB9lF7gAAH0JJREFUnnb+XraXzS01BGu9uSlpKlfGZeCrdg9Bd3aF7DGZ2F5Vz6aKGnZU19NvseKn1ZAbG8nc+GjmxUcTrnPtfupCuAYj/DEcIeQKMhbnH+xs5oXyfexsqydY6831iVP4RlzmpI/Qx7NCmixW9tQ18lVlLTuqG2jqtS/zmhwSxNz4aOYmRDE1Uo/WxaJ1IVyDEf4YjhByBXGG8ws6mnilooDtrXUEar24Nj6Hq+OzCNB6OcnKiWWiKqQsy1S2d7Gjup4d1fUUNrRgsdnw1WrIjYlkbnwUubGRxATqFG8wFcI1GOGP4XiskP923yt8I2wpfmrlGg2d6fyizhZeqShgS0sNOo2Wa+KzuTYhm0Dt5Fq1XKmLtM9kZm9dEzuq69leXU9Tjz1aj9D5kRsTSW5sBLmxkYT6Tfw0CkK4BiP8MRyPFfJLvryfQLUf1+rPI0+XpUjUNR7OP9TVyqsVhWxsrsJPreGq+CyuS8ghZJL0cnGFi1SWZao7u9lT28ie2kb21TXRY7IvvJsUEkhubCSzYiKZER2Oznv8n3xcwSeuhPDHcMbbJ7Is02e2kJ6UCK4k5FuP7OH15s+oGmhkql8K1+vPJ8Ir5PR7OpHxdH55TzuvVhSyrrESL5Wai2JSuS4hh3g/56zSMl644kVqtdkoa+04LuwHGloYsFqRgOTQIKZG6ZkaqWdqVDiROj+nBwWu6BMlEf4YjrN9YrHZKGtp50BjCwcaWjjY2EJSSBD/uvMGcCUhr6urwyrbWN+5mw9aN2HBygXB87kwZD5eqolpMJyIClnZ28mbRw+yuv4IFtnGOeHx3JA4hWlB4YrnfkdiMlykJquVosZW9tc3caChhaKmVoxmCwB6f1+HqNvFPSUsGI1KNabjTQafTCTCH8MZq0+6B0wUN7VyoMEu3CXNrfRbrABE6vyYGqVndmwUty1fBK4m5MfosHTzdst6dvQUodcEcV34+fx/e2ceJMd53uenp+e+r70PLLALYHESIEDwEkiYICmSAg+TUitSSSVHFStOrCqrLJWjSElFrkoqqchlReVQkmlZjmLRZNoURVOiSAo8dJEQSdznAtgF9r5mdu57prvzx8wOCOKgAOzu7Oz2U9X1dc/27Lzzzte/fvv9rpscPfNuxEJWyJl8ludG+vjJ2GkSxQIb3EE+tWIDdzd2IAo3JjRzST1epIqqcj4S59hkmONTYU5Mhqs9YsyiSE/AS2+jn94GP2sbArR5nNc0OKkefTKf6P64lGvxSb5U4mw4xulQhL7QDKdDEUbjKQAMgkBPwFsNRDY0By+aRnrR5cgv12vlVGaQfwrtZbI4wyZ7N58M7qbJPH+LPtSiQmaVIj8fH+DZ4VOMZZO02px8smMdD7X2LIqui0vlIp1OZTgxFaYvFKFvOsLZcKQa4TjNJtY2lIV9TYOfnoCXZpfjik9IS8Unc4Xuj0u5kk8KJYXBaJyzM1FOT0foC0U4H4mjVHQ36LCxtsHPugY/vY0B1jX6sZmurAN1IeQAJU3h9dh+fhp5i5JW4l7vLezx34HVMPe9P2pZIRVN5TehEZ4ZOsmxeAiHaGJPaw9PdPTSbp+7hXKvlaV6kSqqylAsUb2Y+kIRzs3EqheUw2RiVcBLT2XrDnrp8nowG8Ul65PrRffHpQSDQQZGx+gPR+mfiTFQ2YZiCdRKHasGEI2BytOhn6Dj2nph1Y2QzxIrpXh+5pfsSx7HIzp4PLCL21wb53S+jsVSIY/HQ/zz8CnemB5C1TTuCLbz8Y5ebvG3LPj8JIvFJwtBvlTifCRevfD6Z6Kci8Sr+XZREOj0ulnb0kiL3UKn102n1027x7noBi0tJMupjlyOZL7AUDTBSCzBYCzBULRcTifT1XMaHDa6A74LgUHAOydjIupOyGc5lxvn2dBezucnWGFpRgruZo2tY06MWGwVMpTP8MLoGX4yeoZYMUeH3cXj7Wt5sKV7wfqjLzafLDSqpjGRSNE/E6sIfJTheIrxeLJ6jkEQaHU7WeF10+lz0+l10eZ20up24rNZF2Uj9lyyHOpISVWZTKYZT6QYjScZjiUYjiYZiiWIZnPV80yigU6vm3UtjbQ7bfQEfHQHvHis83O91q2QQ/nieid5gudnfkVMSXKzYy1PBHfRaLqx7oqLtUIWVIU3p4b48ehpjsdDWAwi9zWv5NG21ax3z++iyIvVJ7UkGAwyMjHBSKx8IQ/HkgxHEwzHEozGk9X0DIDVaKTV7aDV7aTN7aLV7aDN46LJ5aDBYaurycKuxFKoI5qmkS4UmUplmEimGIunGE9c2KZSmWpKBMBhNpVv3N7yjbvT52aF102zy4FoMCzfAUHXM0Q/rxbZG3uXl6O/Q9EUdnm28jH/nbjE65t0qR4q5OnEDM+Pnmbv5CA5tUS308vDrat5oGXVvETp9eCTheZqPimpKhOJdFUAxhJJxhNpxhNJJhJpiqp60fl+m5VGp51Gp50mp51Gp6N67Ldb8dusiDfYXXK+qYc6ki2WiGRzhNMZplMXtqlUhulUmulUhkwlhTaL22Kmxe2sPl21up20ecrHH/akpQv5dRArJXkx8lt+mzjKKmsrX23/7JU/SNOWRG+EdKnA3slBXhw/S19iBrPBwK7GFTzatpot3qY5i9LryScLxfX6RFFVwplsOcJLzopJuiIm5S2vKBe9RwA8VktZ1O02AnYrPpuVgN2G12bBbbXgsZhxWy24LWZsJuOCp3JqUUcUVSVVKBLP5UnkCiTy5TKSzRHJZIlkchftf1CkoezXi2+idhoddloqou26gVHCupDfAOOFMFk1T7e17YrnDOeneC78Btud67jLs+Wiv9WraJ1JRnhx7CyvTpwjrRTpsLt4uHU1D7V047fc2Fwk9eqT+WS+fKJpGol8galUmlAqy8ysIGWyRLI5ZjIXhEm5wnVqMhhwW824LRbcVjN2kwmH2YTNZLy0NJVLsyhiNorlUhQxGw1YqvsioiDcUPSpaholRaWgKBQUhXxJoXDJsUKmUCJTLJa32f1CkUyxvJ/KF6uCncwXuJJS2U3GCzc+mxWfvXzj81fKWfG2Go3X8vNcE7qQzzNj+RADuTGeCe3l6x2fo93SWP1bvYtWTinxxtQQL46d4Wg8hCgI7PC38kDLKnY2dGAVr73i1rtP5oNa+0TVNJL5ArFsvipsiVyeeL5czu4n8wUyhSLpQpFssUS6WKSoqB/+AZdBoNyoW90MAgbKpWgwUFIUVE1D1aiUF7brxWYyYjeZsJuM2M0mnGYT7tmnD6sZT+UpxP2+0mezYjPNn0D/vtRayGvvgXmmzdLAeCHMGltHVcTzapGDqdPc4d5aY+tuDKto5KHWbh5q7WYwHePn4wP8YvI8/+X4b7CLJnY1dvJAyyq2+poW1ehRnWvDIAh4rJbr6glRVJRyhFuJdLPFUjUyLigqRUWhUFLIz75WUlE0tSzKqoYK5VJTUSpCbbFYKeTzVYEXK2IvCBf2jaIBs2jALIpY3h/9V14zG8WLRNtmMtbtUoCLgSUfkRfVEt+ekNnp3sKtrvUcSPXRlx0mWkwwrcb4iGMT9/tuXRBbFgJV0zgUneLVyXO8OTVEWinSYLFzf/NKPtq8im7n1RdErnX0uRjRfXIxuj8uZUlE5JIkfQX4JtAgy/Ki+oWPZ85RUEvc6loPwAszv+axwF1sa+hl3BTlrYlD1XOv1jBaLxgEgW3+Zrb5m/ny2h38NjzKqxPneHb4JE8PnaDL4WF3Uxe7m7rocnhqba6Ojs4ccMNCLklSB3AfMHzj5sw9byWOcourF4DXY/uxGSxsc5aPGyw+RvMhcmoeAQO/ih9iIDfGRscqdrpvqqXZc4JFNFZFO1bI8cb0EK9PDvKDc0f4+3NH6HH6qn+v5bQAOjo6N8ZcROTfAv4C+Jc5+F9zQlbN8+v4YZpMPlJqlvu8OwD4ZfwQTwR3Vc/7bfgwTtGGAQMvRn5LpJTgI+7NvBR9myaTjzW2zhp9g7nHa7byePtaHm9fSyif4c2pIV6fGuRvBw7xtwOH6HUF2N20giesW6mPJTB0dHRmuSEhlyTpEWBMluUjkiR92LlfAL4AIMsywWDwRj76qmRKOZKpHC+FXsJqsBAIBBjOThKweviDzh0ICBgEA68c3Mc31n2B/bHTOO1OHm7YRZutgYJV5WR2mDuCN8+bjbUkCKxr6+Tfs5PxdIKXR87w0vBpnuw/yJP9B+n1NrC7rZt727pZ72us+3TTjWI0Gue1vtYbuj8updY++VAhlyTpNaD5Mn/6OvA14P7f54NkWX4KeKpyqM13w8DH3XfzMedtvBR5mx+efZFbnOvwCy7Gpiewi1beiB2g3dqIKQXvTB/jYf+dmFIQTofZN3WEDY5VhMNhVE27qDVd1VQMS6gHiBl4NNjFo8EuxjJJ9mdmeGWwj++eeIcnT/yOJoudnQ0d7GzoYKuv+YYXaahH9Ma9i9H9cSkL3Nh5Cdfda0WSpE3A60Cm8lI7MA7skGV58kPevmC9VmZRNY2nJl8gqWRoMvspaQqf6LqPkZlxjqb7ecB3O16jk/FCmGdDr/GZxo9W53KZFfMj6X76MkOMFKZ4PLCLVdbLO7Wema2Q0UKOt8Kj/CY0wrsz4+RVBafRxO2BNu5saOdWfyueOlmH9EbRhetidH9cSt32WpFl+RhQHV0jSdIgsH2x9VqZxSAI/EnLH3IwdZpIKcEdrk10OtoYDI9S0Eo4xfKIyF/HD7HS2oJbdFz03slChGdCv+CzjQ/SbmnghZlf8/mmPXiNzlp9pXnFZ7ayp7WHPa095JQS70Um+PX0MG+Fx9g7NYgBgQ2eILcH27g90MYal3/Zp2B0dGrFkh8Q9EFudq696NhiMBEpJTAKIoO5CY6mB/hS2yexGi7Mu5BRcrwZP8B25zo22FdSVEvsT/VR1C7M53AmO8xqa8eSFDOraKymV1RNoy8xw9vhUd6eGeOpgcM8NXCYoNnGbRVR3xFowWGc/9XtdXR0ysyZkMuy3DVX/2sh6bA00Wzy85Xzf8NmRw/3+3ZUUyqz/cqnilGmi1E+2/gAAFElid/oRtEUhnKT7E/1cTY3QlrJ8aDvNu5wb6rlV5pXDILAek+Q9Z4g/6Z7C5F8lt/NjLNvZpRfTg/xs/F+REFgvTvILf4Wtvlb2OgJYjLU//StOjqLlWUXkX8Qm8HCv2q4j4f9OxEFAYtQjiTH8iHaLA0AHM8M4De6K+KtEirGiJWSNJsDfHP0aW52ruWr7Z/lQKqPgdxYLb/OguO32KrTBJRUlWPxEO/MjPNeZJz/c/4YPzh/FKvByE2+Rrb7mtnub2G1y68Px9bRmUOWvZDP4hAvNNwlSmmeCe/lK22fBuBE5jyP+ncCEC0lOJYe4GbnWs5kR4gpKXZ7t6NpGp2WZvYlj5NWsjhEG+8lT3E+P85WxxpWz9GqRosZo8HAVl8TW31N/AlbSRYLHIxOciAywf7oJE/2HwTAbTKz1dvMFl8jW7xN9Lh8+lwwOjo3gC7kl8FtdFRFXNM0tjrWcC43To+1nddi+ylpCre5NvKt8Wf5qLc8T4sgCJzKDKJqGnaDlZ9H9nE2N0KvbQU/nH6ZP2p8iB5bey2/1oLjMpm5u7GTuxvLA6tCuQz7oxMciExyODbFr0LlwcAO0cRmbyM3+RrZ6m2i1x3QUzE6OtdAzYQ8r8WwCN5affzvjSAIbLSv4kehVzma7me1rYOHfLcTV1Lk1AJ3uDdVc+mvxt7hc40PsS95nISS5jH/XaywNuMUbRxN9y87If8gDVY7D7Z082BLNwDTuTSHY9MciU5xODbFvv5yWspiENngCbLR08gmTwMbPEG8y6Sro47O9VAzIX8r80U6TXvoMj2CUbi+pdoWijZLA/+h/TPESylcoh2DYGCmGKfD0khJUzAaRN5LnsJqMLPG1sHL0X3c672lmmM/kTlPm7m8/8EBRsuZRquD+5tXcn/zSgCihRxHY9Mcjk5xJDbN00PHqwsqdNhdbPQ0sNHTwCZPIyudHj0do6NToWZCHhRv5nzxOUaLr9BleowO04OIwsKsFH+9eIxOZgdQ+Ywu4qU0/zj9Cs1mP+OFME8E/oCz2RF8Rhet5iBGQSRaSpJQ0uxx3gmgi/hV8JmtF6VickqJU4kZjsdDHI+H2Bce5+WJcwDYRCNrXQHWuvysdZfLTodbF3edJUOmVORsKsLpRITTyRm8Jiv/vfWRy55bMyHfbP1zEsqj9Bef4WzxRwyVfsZK0+O0G+/DIJhqZdaHMttP3CAY+Hctf8hrsf0klTSPB3bRYPJyNN2PKIhYKt/hN4kjtJsb8YpLc+DQfGIVjdXGUyi3V4xlUxyPhzgZD9GXjPDC2BnyI+V1L22ikdVOH2vcAXpdAda6/aywe5bltAI69UWqVOBssizafckZTiciDGfi1eXt/GYrdwavnJqtaWOnW+zmZvE/EVVOMVB4htOFHzBUfJFVpk/QYrwbg7C422JNgpEHfbdd9JrFYGamGMcuWomWkryTPMEfNz2CXbx8jjenFi4afKRzZQRBoN3uot3u4oGWVUB5FfuhTJwzlQvgTDLCS+P9PKf0AWAUDKx0eOh2+uh2eul2+ehx+giYbdc8eGvv5HnenB6i2erk1kArtwbKUzTECjk8JsuSHAymM7eUVJWRTIL+VJSBVIyBVJT+VJSpXLp6ToPFzlqXn3ubu+h1BVjj9tNguXr6eVEopU9cxzbrXxJRj9Ff+CdOFr7L+eKPWWl6oi4E/f00m/xYDCa+Ovgd1tg62e7spcvactlzQ8UY3xj+Plsca7jDvYl1thVLakKuhcBoMFRE2seDlBtRFa18sZxORBioXDAHopO8Mnmu+j6PyUK308cqp5cuh4eVDg8bPQ1X7C1zLBZCHjnFH63cTF9ihp+Nn2WV08vpxAzf6T+IqmnsbOjgT1dvu+z780oJk0HUU2vLhJKqMp5NMZiOMZiJM5SOM5CKMZiOUVDL66iKgsAKu4ebvI10O8sBxlqX/7oWUF80CikIAgFxM37rJsLKAQaKcl0Kusfo5N82P8ZIfgqnaMdTSalcrpFTxMCd7s28kzzJu6mTeEQnt7k2cJtrw0WLROtcG6JgoMvhpctxca+oRDFfjoCS0WpE9NJ4P1mlPNXCa7s+dUUhH8rE0TS4M9jONl8z/+3kW7w8McDR2DR/teUerAYj3z7zHifjYdZ7Lp7ONJzP8L3+Q7wbGafL4eE/r/8IDdZyhDWVS3MoOkWn3c1ql0/vdllnpEsFRjNJfpcOc2xylMF0nMF0nNFMgqJ2YeHroMVGt9PHLf51rHJ66XH6WOHwYJ6j33vRrtmpaRph5SDnijIJdQCb0EiX6XFajXfPWQ59scziVtRKHE33sy95nOPpcyiodJgbud29kVuc6xd0Yq7F4pOFQtM0pvMZRjNJtvkvN1tz2SdHhs/znf6DzOQz+Mw2tviaAI0DkUn+55Z7yCsl/t/IKUqqyudX3YSiqYiCgYKq8NxIH2eSEb6xcSfy8ClOJsJ8Y+NOBtNx/uHcEdJKkWypxH3NK3msfQ0A/ckoPzh3hIDFxn3NK9nsXTw39uVURzRNI17MM5ZNlrdMitFsgtFMktFskmghVz3XgECrzUmXw8MKh4euyrbC4cE5R3MPzeuanfOBIAg0GLcRFG+uCvqpwvc4X3yOLtOjtBrvWfS9XH5fTIKRbc5etjl7SSoZ3kue4u3kMeTwG/xz+E3W2VZwq2sDWxyrr5hr17k+BEGgyeqgyeq46nn7ZsYIWGzc07SCX00P4zVZmMql6XS4gfKyegVFoaCWG141DRDgbDJCKJdhT2sPUM5/zgZPb4dHUYG/2rKbQ9FJ/nHwOI+1r6EvMcPzo6e5yddEopjn5YkBet2BOYvedC5QUlVC+QxTuTSTuXSlTDFV3U9Xn9hmabDY6bC7+EiwnXa7m3abi81tnTjzKhaxNr/RohXyWd4v6DPKYc4Vn6Ov8PecK/6YTuPHaDfdj0m4+kVYT7hEO/d4t3GPdxsThRneTZ7gneRJ/mH6JYyIbLCvZLtrHTc5erAZlsaNbLEznU1xPhXjnqYutvmbcRhNHI+V55TZ09ZTPS+Uz9DjrEy4VulvMJ3LVG8WUM7fN1odnIiHSBTz1QZTh9FMh93NcDrBkdg0eaXEJzvXkSoV+O7Zg+wLj3F3Y+eSWCB8IcgrJaLFHNFCjnA+y0w+SzifIZTPMFPIEq4cRws5PpiT8JqsNFsddNo97PC30mR10GZ30WZz0WZzYhEvlc2gt7ZPKYteyGcRBIGgcSsBcQtR9QSDxZ/QX3ya88Xn6TDdT6dxDxaDr9Zmzikt5gCPBu7iEf9OBvMTvJc6xf5UH0em+jEKIpvs3Wx39rLZ0aP3fJlHrKKR/lSUR9pWA5BXFBRNI5TP4DWVn5Ai+SwZpUinwwNQFYd0qYBVFLFUoul4MY9dNJFVSuRVheaKwAuAomlklCKT2RSrXX4A7KIJgyAQzme4HEVVoaCq2EXjkhR4RVNJl4okigVSpQLJYp5kqUiylCdZLBAv5okWcsSKOWKFfKXMkflAFA1lH/vNNgIWG0GLjV53gKDZRqPVQXPlqazJ6sB6GaFe7NSdxYIg4Bc34hc3klDOMVh8gcHiTxkqvkSrcRcrTI/iMFy+l0i9IggCK62trLS28vHAPZzLjbE/dYr9qdMcSp/BLBjZZO9mm7OXTY5VWPVIfU5xm63cFmjjb87up83mYjyb4uMdvXjNFk4mwmz3t/DziQEaLA5WVyLy2TRIs83JiUS4KuxDmQRdDg9uk4VUsYCrkjtNl4oomorLaCaUz7CpkhM3CAKhfKbal/6DHI5O8WeHXkMUBJxGM3bRhFU0YheN2IxGrAYjdqMJm2jEKpaPzQYRk8GAySBifn8piJgNIqIgIAhCuUTAIJTtMFB+3ScUicZjqJqGqmloaCiahqZpqJQb9kuaWr3JlEuFoqpS0CqlqpBVSuSUElmlSKY0u18ioxTJKSXSpSJppXjV38ZsMOA1WfGarXhNFtrtLnyzx2YLPpOVoMVOg8WOz2xdsmMK6k7I349bXMVm8c/JqJMMFX/KeOkNxkqv0yjeSpfpMTxiz4f/kzrDIAj02NrpsbUjBXfTnxvlvVQfB1J9HEifxiiIbLCtZKtzDTc5VldXPtK5MT61Yj3b/M2MZ1PcHhSraY5vHP8Nn3/3Z7RYnXx1/e34LTZOJ2YYyyb5SEMHW7xN/N3AYSKFLE1WB+/OjPNA8yrWuPz0JWewieWG+xPxEC6jmTa7i5FMgqZKr5asUiSrlGi2lhu8Pxh1t9ldfHH1tmrEmq00nM6KYbyQJ1MRy1mhrEn3hg8gCgJ28cINZnbfb7HSJjqxiSbsogmXyYzTaMZtNOM0mXEZzdXXXEYztiX6JHKt1LWQz2I3NLPO8sesMn+CkeLPGSm+woxymLvsf4dRWLpCZhAMrLF1ssbWyaeC99KfG+Vg6gyH0mc4Mt3PSksLX+v4XK3NXBJYRSObvY0X9R4RBIG/3HQXUG40m4323CYLM4UsJsGAIAg80dHLfz3xFgVV5TNdG+h1BwDY7m/hp+Nn2epr4vXpIb68dgcAbTYX4XwWgLdCo/jNNtpsl++51Gpz8ekVG37v76Fp5ei5HCErFCqRc1FVq5GyinbVaNvpdpFOJK8atRsFAyaDoRr9mwXxoicAvT/93FKz7oe1+FAdHR2dJcAld8FaJYyExbBJknSg1jYstk33ie4T3R+L3ieXsDQz/zo6OjrLCF3IdXR0dOqc5S7kT9XagEWI7pNL0X1yMbo/LqWmPqlVY6eOjo6Ozhyx3CNyHR0dnbpHF3IdHR2dOmdJDAiaCyRJ+grwTaBBluXlMUfnFZAk6ZvAw0ABGAD+tSzLsdpatfBIkvQA8G1ABL4vy/L/qLFJNUWSpA7g/wLNgAo8Jcvyt2trVe2RJEkE9gNjsizvqYUNekROtYLeBwzX2pZFwl5goyzLm4EzwH+ssT0LTuXifBJ4EFgPfEqSpPW1tarmlIAvy7K8DrgN+FPdJwD8GXCqlgboQl7mW8BfoI84BUCW5V/Isjw7fdzvgCuv+rp02QH0y7J8TpblAvAs8GiNbaopsixPyLJ8sLKfpCxebbW1qrZIktQOfAz4fi3tWPZCLknSI5QfiY7U2pZFyueBl2ttRA1oA0bedzzKMhet9yNJUhewFXinxqbUmv9FOQhUP+zE+WRZ5MglSXqNcl7vg3wd+Bpw/8JaVHuu5hNZlv+lcs7XKT9OP72Qti0SLjcUWn9iAyRJcgI/Br4ky3Ki1vbUCkmS9gDTsiwfkCRpVy1tWRZCLsvyvZd7XZKkTcBK4IgkSVBOIRyUJGmHLMuTC2jignMln8wiSdLngD3AblmWl6OAjQId7ztuB66+0OwyQJIkE2URf1qW5edrbU+NuRN4RJKkhwAr4JYk6UeyLH9moQ3RBwS9D0mSBoHteq8V6QHgr4G7ZVkO1dqeWiBJkpFyQ+9uYAx4D/i0LMsnampYDZEkSQB+CERkWf5Sre1ZTFQi8q/ovVZ0FhP/G3ABeyVJOixJ0vdqbdBCU2ns/SLwKuVGPXk5i3iFO4HPAvdU6sXhSjSqU2P0iFxHR0enztEjch0dHZ06RxdyHR0dnTpHF3IdHR2dOkcXch0dHZ06RxdyHR0dnTpHF3IdHR2dOkcXch0dHZ065/8DwfS/30q2ijIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "contour = ax.contour(x_mesh, y_mesh, z_mesh)\n",
    "ax.clabel(contour, inline=1, fontsize=10)\n",
    "ax.set_title('Линии уровня функции fun')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Задачутка 2. Велика ты, правда, широка стоишь! Ты горами поднялась до поднебесья!\n",
    "\n",
    "(с) Алексей Толстой\n",
    "\n",
    "[Как известно](https://www.livelib.ru/quote/305456-malysh-i-karlson-kotoryj-zhivet-na-kryshe-astrid-lindgren), Фрекен Бок пьёт коньяк по утрам.\n",
    "\n",
    "![](http://semyarf.com/UPLOAD/2016/01/21/frekenbok-216_700_0.jpg)\n",
    "\n",
    "А у нас даже есть дневные данные (в граммах):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = [3.2, 7.9, 5.4, 4.9, 6.2, 4.3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "а) Предполагая, что $y_i$ независимы и нормальны $N(\\mu, \\sigma^2)$ выпишите функцию правдоподобия.\n",
    "\n",
    "Чтобы гарантировать положительность параметра $\\sigma^2$ функция у нас будет зависеть от вектора параметров $\\theta$, причем $\\theta_1 = \\mu$, а $\\theta_2 = \\ln \\sigma^2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def minus_lklh(theta, y):\n",
    "    # theta_1 = mu\n",
    "    mu = theta[0] \n",
    "    # theta_2 = log(sigma^2) \n",
    "    s2 = np.exp(theta[1]) # лайфхак для неотрицательности sigma^2\n",
    "    n = len(y)\n",
    "    lklh = - 0.5 * n * np.log(s2) - 0.5 * np.sum((y - mu) ** 2) / s2\n",
    "    return -lklh\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "б) Найдите экстремум функции правдоподобия"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.5819683626297953, 3.5000002497146308, 0.2750002585914815)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "theta_init = [0, 0]\n",
    "y = [3.4, 2.7, 3.8, 4.1]\n",
    "\n",
    "optim_res = minimize(minus_lklh, theta_init, args=y)\n",
    "optim_res.fun, optim_res.x[0], np.exp(optim_res.x[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "в) Получите тот же результат используя встроенный метод `stats.norm.fit`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3.4999999999999996, 0.2749999999999998)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu, sigma = stats.norm.fit(y)\n",
    "sigma2 = sigma ** 2\n",
    "mu, sigma2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Задачулечка 3. \n",
    "\n",
    "А ещё Фрекен-Бок иногда видит привидения! Данные по количеству привидений у нас тоже есть :)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "fbock = pd.DataFrame({'y': [3.2, 7.9, 5.4, 4.9, 6.2, 4.3], 'ghost': [1, 2, 0, 0, 2, 0]})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Предположим, что количество привидений имеет пуассоновское распределение с параметром $\\lambda$. \n",
    "\n",
    "а) Оцените $\\lambda$ с помощью ММП:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5.9116077839764465, array([0.83333469]))"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def minus_lklh(theta, ghost):\n",
    "    rate = np.exp(theta) # тот же трюк, чтобы лямбда никогда не стала отрицательной\n",
    "    log_probs = -rate + ghost * np.log(rate) # про минус log(k!) забиваем!\n",
    "    lklh = np.sum(log_probs)\n",
    "    return -lklh\n",
    "\n",
    "theta_init = 0\n",
    "\n",
    "optim_res = minimize(minus_lklh, theta_init, args=fbock['ghost'])\n",
    "optim_res.fun, np.exp(optim_res.x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Предположим, что в $i$-й день интенсивность пуассоновского распределения $\\lambda_i$ связана с количеством выпитого коньяка формулой $\\lambda_i = \\exp(a + b y_i)$.\n",
    "\n",
    "б) Оцените параметры $a$ и $b$, лапками выписав функцию правдоподобия. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4.901374802427522, array([-2.59878425,  0.41696008]))"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def minus_lklh(theta, fbock):\n",
    "    a, b = theta # распакуем параметры для читабельности\n",
    "    rate = np.exp(a + b * fbock['y']) # тот же трюк, чтобы лямбда никогда не стала отрицательной\n",
    "    log_probs = -rate + fbock['ghost'] * np.log(rate) # про минус log(k!) забиваем!\n",
    "    lklh = np.sum(log_probs)\n",
    "    return -lklh\n",
    "\n",
    "theta_init = [0, 0]\n",
    "\n",
    "optim_res = minimize(minus_lklh, theta_init, args=fbock)\n",
    "optim_res.fun, optim_res.x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Данная модель является довольно популярной и реализована в `statsmodels`.\n",
    "\n",
    "в) Сравните ваши результаты с результатами в `statsmodels`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "# !pip install statsmodels\n",
    "import statsmodels.formula.api as smf\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 1.047945\n",
      "         Iterations 6\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Poisson Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>         <td>ghost</td>      <th>  No. Observations:  </th>  <td>     6</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                <td>Poisson</td>     <th>  Df Residuals:      </th>  <td>     4</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>                 <td>MLE</td>       <th>  Df Model:          </th>  <td>     1</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>            <td>Mon, 25 May 2020</td> <th>  Pseudo R-squ.:     </th>  <td>0.1384</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                <td>01:38:59</td>     <th>  Log-Likelihood:    </th> <td> -6.2877</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>converged:</th>             <td>True</td>       <th>  LL-Null:           </th> <td> -7.2979</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>     <td>nonrobust</td>    <th>  LLR p-value:       </th>  <td>0.1552</td> \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>   -2.5988</td> <td>    1.918</td> <td>   -1.355</td> <td> 0.175</td> <td>   -6.358</td> <td>    1.161</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>y</th>         <td>    0.4170</td> <td>    0.297</td> <td>    1.404</td> <td> 0.160</td> <td>   -0.165</td> <td>    0.999</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                          Poisson Regression Results                          \n",
       "==============================================================================\n",
       "Dep. Variable:                  ghost   No. Observations:                    6\n",
       "Model:                        Poisson   Df Residuals:                        4\n",
       "Method:                           MLE   Df Model:                            1\n",
       "Date:                Mon, 25 May 2020   Pseudo R-squ.:                  0.1384\n",
       "Time:                        01:38:59   Log-Likelihood:                -6.2877\n",
       "converged:                       True   LL-Null:                       -7.2979\n",
       "Covariance Type:            nonrobust   LLR p-value:                    0.1552\n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept     -2.5988      1.918     -1.355      0.175      -6.358       1.161\n",
       "y              0.4170      0.297      1.404      0.160      -0.165       0.999\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "smf.poisson(data=fbock, formula='ghost ~ 1 + y').fit().summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Задачушка 4\n",
    "\n",
    "Предположим другую вероятностную модель. Теперь мы будем предполагать, что сам факт обнаружения хотя бы одного приведения, имеет вероятность $p_i$, зависящую от $y_i$:\n",
    "\n",
    "$p_i = \\exp(a + b y_i) / (1 + \\exp(a + b y_i))$\n",
    "\n",
    "а) Введите дамми-переменную равную 1, если Фрекен Бок видела приведение:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "fbock['atleast1'] = 1 * (fbock['ghost'] > 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "б) Определив лапками функцию правдоподобия, оцените $a$ и $b$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3.8709632618137633, array([-2.35157624,  0.44383861]))"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def minus_lklh(theta, fbock):\n",
    "    a, b = theta # распакуем параметры для читабельности\n",
    "    p = np.exp(a + b * fbock['y']) / (1 + np.exp(a + b * fbock['y']))\n",
    "    log_probs = fbock['atleast1'] * np.log(p) + (1 - fbock['atleast1']) * np.log(1 - p)\n",
    "    lklh = np.sum(log_probs)\n",
    "    return -lklh\n",
    "\n",
    "theta_init = [0, 0]\n",
    "\n",
    "optim_res = minimize(minus_lklh, theta_init, args=fbock)\n",
    "optim_res.fun, optim_res.x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "в) Оцените $a$ и $b$ с помощью `statsmodels`. Нужная модель называется `logit`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.645161\n",
      "         Iterations 5\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Logit Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>       <td>atleast1</td>     <th>  No. Observations:  </th>  <td>     6</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                 <td>Logit</td>      <th>  Df Residuals:      </th>  <td>     4</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>                 <td>MLE</td>       <th>  Df Model:          </th>  <td>     1</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>            <td>Mon, 25 May 2020</td> <th>  Pseudo R-squ.:     </th>  <td>0.06923</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                <td>01:45:29</td>     <th>  Log-Likelihood:    </th> <td> -3.8710</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>converged:</th>             <td>True</td>       <th>  LL-Null:           </th> <td> -4.1589</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>     <td>nonrobust</td>    <th>  LLR p-value:       </th>  <td>0.4479</td> \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>   -2.3516</td> <td>    3.376</td> <td>   -0.697</td> <td> 0.486</td> <td>   -8.969</td> <td>    4.266</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>y</th>         <td>    0.4438</td> <td>    0.620</td> <td>    0.716</td> <td> 0.474</td> <td>   -0.772</td> <td>    1.660</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                           Logit Regression Results                           \n",
       "==============================================================================\n",
       "Dep. Variable:               atleast1   No. Observations:                    6\n",
       "Model:                          Logit   Df Residuals:                        4\n",
       "Method:                           MLE   Df Model:                            1\n",
       "Date:                Mon, 25 May 2020   Pseudo R-squ.:                 0.06923\n",
       "Time:                        01:45:29   Log-Likelihood:                -3.8710\n",
       "converged:                       True   LL-Null:                       -4.1589\n",
       "Covariance Type:            nonrobust   LLR p-value:                    0.4479\n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept     -2.3516      3.376     -0.697      0.486      -8.969       4.266\n",
       "y              0.4438      0.620      0.716      0.474      -0.772       1.660\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "smf.logit(data=fbock, formula='atleast1 ~ 1 + y').fit().summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}