Elistratov 382008-2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 29. Нелинейный метод наименьших квадратов\n",
    "Нелинейным методом наименьших квадратов\n",
    "построить модели вида\n",
    "$$\n",
    "y = \\frac{\\beta_1}{1 + \\beta_2 e^{-\\beta_3 t}},\n",
    "\\qquad\n",
    "\\mbox{и}\n",
    "\\qquad\n",
    "y = \\frac{\\beta_1}{(1 + \\beta_2 e^{-\\beta_3 t})^{1/\\beta_4}},\n",
    "$$\n",
    "для данных, приведенных в таблице\n",
    "```\n",
    "    t         y    \n",
    "  ----------------   \n",
    "    1       16.08    \n",
    "    2       33.83  \n",
    "    3       65.80  \n",
    "    4       97.20  \n",
    "    5      191.55  \n",
    "    6      326.20  \n",
    "    7      386.87  \n",
    "    8      520.53  \n",
    "    9      590.03  \n",
    "   10      651.92  \n",
    "   11      724.93  \n",
    "   12      699.56  \n",
    "   13      689.96  \n",
    "   14      637.56  \n",
    "   15      717.41             \n",
    "```\n",
    "Предиктор $t$ – время, \n",
    "переменная отклика $y$ – вес сухой луковицы вместе с надземной частью растения.\n",
    "Найти общую ошибку в каждом из этих случаев и сравнить результаты."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math as m\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import curve_fit "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Исходные данные\n",
    "t = np.arange(1, 16, 1)\n",
    "y = np.array([16.08,33.83,65.80,97.20,191.55,326.20,386.87,520.53,590.03,651.92,724.93,699.56,689.96,637.56,717.41])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Метод наименьших квадратов (МНК).  \n",
    "Пусть имеются точки\n",
    "$$\n",
    "(x_1, y_1),\\, (x_2, y_2),\\, \\dots,\\, (x_m, y_m).\n",
    "$$\n",
    "Ищем линейную функцию  \n",
    "$$\n",
    "y = \\beta_0 + \\beta_1 x,\n",
    "$$\n",
    "такую, что\n",
    "$$\n",
    "y_i \\approx \\beta_0 + \\beta_1 x_i\n",
    "\\qquad\n",
    "(i=1,2,\\dots, m).\n",
    "$$\n",
    "А именно, ищем $\\beta_0$, $\\beta_1$, на которых достигается минимум суммы квадратов невязок (*остаточной суммы квадратов*)\n",
    "$$\n",
    "R(\\beta_0,\\, \\beta_1) = \\sum_{i=1}^m \\bigl(y_i - \\beta_0 - \\beta_1 x_i\\bigr)^2 \\to \\min_{\\beta_0,\\, \\beta_1} \n",
    "$$\n",
    "$R(\\beta_0,\\, \\beta_1)$ - это сумма квадратов отклонений предсказанных значений $\\beta_0 + \\beta_1 x_i$ от истинных $y_i$. \n",
    "\n",
    "Для построения модели с помощью метода наименьших квадратов используем функцию curve_fit, которая использует нелинейный метод наименьших квадратов для максимального совпадения построенной функции с исходными данными.\n",
    "\n",
    "Построим модель вида $ y = \\frac{\\beta_1}{1 + \\beta_2 e^{-\\beta_3 t}} $, применив МНК.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[β1 β2 β3] = [7.02871478e+02 8.49924183e+01 6.88565587e-01]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXhU1fnA8e+bhAQSwp6wZBsQUGTHgFitda8iirXYHxIVlRqtUveFGncbi3ZBrUubuqFE3FFUakUUrbUuYQtL2MQkhAAJW8CEJcv7+2NuaIBJCDDJnZm8n+eZZ+4998zkHci8Offcc88RVcUYY0xoCXM7AGOMMf5nyd0YY0KQJXdjjAlBltyNMSYEWXI3xpgQFOF2AABdunRRj8fjdhjGGBNU5s+fv1lV43wdC4jk7vF4yMnJcTsMY4wJKiJSUN8x65YxxpgQZMndGGNCkCV3Y4wJQQHR5+5LZWUlRUVF7N692+1Qmlzr1q1JTEykVatWbodijAkRAZvci4qKiI2NxePxICJuh9NkVJUtW7ZQVFREz5493Q7HGBMiArZbZvfu3XTu3DmkEzuAiNC5c+cWcYZigk92djYej4ewsDA8Hg/Z2dluh2QaKWBb7kDIJ/ZaLeVzmuCSnZ1Neno6FRUVABQUFJCeng5AWlqam6GZRgjYlrsxxl0ZGRn7EnutiooKMjIyXIootDT1WZEl9waICJdffvm+/aqqKuLi4hg9ejQAL730EmFhYeTm5u6rM2DAAPLz8wHvzVmbN28GYP78+fTs2ZOFCxc23wcw5igUFhYeVrlpvNqzooKCAlR131mRPxO8JfcGxMTEsHTpUnbt2gXAnDlzSEhI2K9OYmIimZmZDb5Pbm4uY8eO5fXXX2fo0KFNFq8x/pScnHxY5abxmuOsyJL7IZx33nl8+OGHAMyYMYNLL710v+OjR49m2bJlrFy50ufr8/LyuOiii3jllVcYMWJEk8drjL9kZmYSHR29X1l0dPQhGzPm0JrjrCigL6jWevD9ZSwv3uHX9zy+Rzvuv6D/IeuNGzeOhx56iNGjR5Obm8vVV1/Nv//9733Hw8LCuPPOO3nkkUeYNm3aQa8fM2YM06dP55RTTvFr/MY0tdqLphkZGRQWFpKcnExmZqZdTPWD5ORkCgoOnhbGn2dFh2y5i8ixIrKozmOHiNwsIp1EZI6IrHaeOzr1RUSeFJE1IpIrIsP8Fq0LBg0aRH5+PjNmzGDUqFE+64wfP56vv/6aH3744aBjZ511Fs899xzV1dVNHaoxfpeWlkZ+fj41NTXk5+f7JbHb8MrmOSs6ZMtdVVcCQwBEJBxYD8wEJgNzVXWKiEx29u8CzgP6OI8TgWed5yPWmBZ2U7rwwgu5/fbbmTdvHlu2bDnoeEREBLfddhuPPvroQceeeuoprrvuOq6//nr+/ve/N0e4xgQsG17p1RxnRYfbLXMm8L2qFojIGOA0p3waMA9vch8DvKyqCnwtIh1EpLuqbvBTzM3u6quvpn379gwcOJB58+b5rHPllVfy2GOPsXPnzv3Kw8LCmDFjBj//+c+57777eOihh5ohYmMCU0MXEgMluVdW11Cxp5of91ZRsaeK8r3VlO+pYndlNZXVNeytViqrathbXePdr92uUue4t6yyznNlte5X7j2WQr9bXuHp03tz/qDufv8ch5vcxwEznO2utQlbVTeISLxTngCsq/OaIqdsv+QuIulAOgT+1ffExERuuummButERkZy4403+qwXFRXFe++9x89+9jO6du3KDTfc0FShGhPQmupCoqry/EuvcP9997BhfRFx3RO4fNJdDD/rQsr3eJNzbZIu31tF+Z4qKmr391TvKyvfW83eqpojjiM8TGgVLrQKDyMqIoxW4d5HpLMdGS77tqMjI2gVHkZMVPhRffb6iLeB3YiKIpFAMdBfVTeJyHZV7VDn+DZV7SgiHwJ/UNUvnfK5wJ2qOr++905NTdUDF+vIy8ujX79+h/+JglRL+7ymZfJ4PD4vJKakpOy7P6SuvVU1lP64h9Kd3kfJzt11tv9X/v1/P6Jk9pNo1Z59r5WIKDqdO4m2/U8HINJJpNGREbSNiiA6KpyYyAhi9j17y9pGRhAdFUFMZDgxURH7XtOmVbiTqIXI8HBaRci+xB3pJPHwsOa921xE5qtqqq9jh9NyPw9YoKqbnP1Ntd0tItIdKHHKi4CkOq9LxPtHwRjTwmVmZu7X5w4Q1boNZ11+E4/MzqNkx+59ybxk5x62V1T6fJ9OMZHEx0YRFxtFry4xLJ86fb/EDqBVe2i96A0Wvf4I0ZERREa0rJHfh5PcL+V/XTIAs4AJwBTn+b065ZNE5DW8F1LLgrm/3RhzdLaW7yVvww7yNuxgeav+9L74NpbP+htVOzYT3q4LbU+9gk8q+xL5VT7xsVHEx0bRs0sMI3p2Ij62NXFOWZzz6NI2ilbh+yfqxy/d6PNnF68vokN0ZHN8zIDTqOQuItHA2cC1dYqnAG+IyESgELjEKZ8NjALWABXAVX6L1hgTsKprlB82/8jyDTv3JfO8DTvYtON/Leq42CgGnz6aceMv5bhusXRr18abvNtFERsVccST6DXHuPFg06jkrqoVQOcDyrbgHT1zYF0F7IqhMSFsx+5KVjhJfHnxDvI27mDlxp3scS5GRoQJvePb8pNjutCveyz9urejX/d2dGkb1STx+Oruael30wbFHarGGPes376LJUXb92uRF23bte94x+hW9OvejstGpjhJPJbe8W2JimiaUSC+2N20B7Pkbow5SPH2XXyYu4H3c4vJLSoDIEygZ5cYhiR14NIRyRzvtMa7tosKiDUJ0tLSWnQyP5Al9waEh4czcOBAKisriYiIYMKECdx8882EhdV/1T0/P5+vvvqK8ePHN2Okxhy90p17+OfSDby/uJjv8rcBMCixPXePOo4Te3amb9dY2kQ2X2vcHJ2QSe7Z2dl+PyVr06YNixYtAqCkpITx48dTVlbGgw8+WO9r8vPzefXVVy25m6BQVlHJR8s28P7iDXz1/WZqFI7tGsvt5/Rl9KAeeLrEuB2iOUIhkdybY76K+Ph4srKyGD58OA888AAFBQVcfvnllJeXA945ZH7yk58wefJk8vLyGDJkCBMmTOAXv/iFz3rGuOXHPVV8snwT7y8u5ovVpVRWK57O0dxwem9GD+rBsd1i3Q7R+IOquv444YQT9EDLly8/qKw+KSkpChz0SElJafR7+BITE3NQWYcOHXTjxo1aXl6uu3btUlXVVatWae1n+Oyzz/T888/fV7++egc6nM9rzIGmT5+uKSkpKiKakpKi06dP3+/4rr1VOju3WH8zPUf7ZszWlLs+0JMe+UQzP1yuueu2a01NjUuRm6MB5Gg9eTUkWu7NuRyYOtM1VFZWMmnSJBYtWkR4eDirVq3yWb+x9Yw5UvWduVbVKEnDz+b9xRv4eNlGyvdW06VtJOOGJ3HB4B4MS+5IWDPfLm+aT0gk9+a6gWHt2rWEh4cTHx/Pgw8+SNeuXVm8eDE1NTW0bt3a52umTp3aqHrGHKn6Zlq85re30eO6F2jfphUXDO7BBYN7cGLPTkSEt6zb8FuqkEjuzXEDQ2lpKddddx2TJk1CRCgrKyMxMZGwsDCmTZu2bzGO2NjY/ab9ra+eMf5S3xlq5Y5SXrgylVN6x7W4eVVMiKyhmpaWRlZWFikpKYgIKSkpZGVlHfXF1F27djFkyBD69+/PWWedxTnnnMP9998PwPXXX8+0adMYOXIkq1atIibGO6pg0KBBREREMHjwYKZOnVpvPWP8QVWJ69bD57GU5GTOOK6rJfajFLQrR9XXGd+cj6O9oBoKWtrnNUcvf/OPOuGFb7Tz6Ns0rFXUfoMJoqOjD7qoag7f9OnTNTo6OmD/bWnggqr9STcmyOyurObxT1Zx9tQvyMnfxmN33cBLzz/n9zNX0/DKUYEuJPrcjWkpPl9Vyn3vLaVgSwUXDO7BPef3o2u71nBKTy6//DK3wws5zTkSz98COrmrakDMWdHUtJGrYZmWa0PZLh7+YDmzl2ykV5cYpk88kVP6dHE7rJAXzFMJB2y3TOvWrdmyZUvIJz5VZcuWLTZE0vhUWV3DP75Yy5l//py5eSXcfk5f/nnzTy2xN5PMzEyio6P3KwuWqYQDtuWemJhIUVERpaWlbofS5Fq3bk1iYqLbYZgA8+0PW7n33aWs3LSTM4+L54EL+5PUKfrQLzR+E8xTCTd6geym5GuBbGNaqs0/7uEPs1fw9oIiEjq04YEL+3P28V3dDssEIH8tkG2MaULVNcqMbwt57KMV7Kqs5vrTjmHSGb2JjrSvqTl89ltjTABYUlTGPe8uYXFRGSf16szDF/Wnd7zNzmiOXGMXyO4APAcMwDuQ/2pgJfA64AHygV+p6jbxDm95Au8i2RXAlaq6wO+RGxMCynZV8qd/rWT6NwV0joniiXFDuHBwjxYxSsw0rca23J8APlLVsSISCUQDdwNzVXWKiEwGJgN3AecBfZzHicCzzrMxxqGqzFy4nkdm57G1fC8TTvJw6zl9ade6lduhmRBxyOQuIu2AU4ErAVR1L7BXRMYApznVpgHz8Cb3McDLzq2xX4tIBxHprqob/B69MUGoYEs5d7yVy7c/bGVwUgdeumoEAxLaux2WCTGNGefeCygFXhSRhSLynIjEAF1rE7bzHO/UTwDW1Xl9kVO2HxFJF5EcEclpCcMdjQH48zPP0/eYXrx53cmUT7uWS9rnW2I3TaIxyT0CGAY8q6pDgXK8XTD18dVZeNB4S1XNUtVUVU2Ni4trVLDGBLM/PfM8d9x8A3vLSgBl88b1XHttevDMMmiCSmOSexFQpKrfOPtv4U32m0SkO4DzXFKnflKd1ycCxf4J15jgtKy4jLvvvhut3LNfebBMQmWCzyGTu6puBNaJyLFO0ZnAcmAWMMEpmwC852zPAq4Qr5FAmfW3m5ZsSVEZ4//xDZVlvrsfg2ESKhN8Gjta5rdAtjNSZi1wFd4/DG+IyESgELjEqTsb7zDINXiHQl7l14iNCSKL1m3n8ue/oV3rViQkJrK+aN1BdYJhEioTfBqV3FV1EeDrFtczfdRV4IajjMuYoDe/YBtXvvAtHWMiefWaE/k86Q9NvhykMbUCdlZIY4LZtz9s5Yrnv6FLbBSvXzuSxI7RTbYcpDG+2MRhxvjZV99vZuJLOfTo0JpXrxnpXUzDmCZgE4cZ00y+XL2ZX7/8HUkdo3n1mpHExUa5HZJpoaxbxhg/mbeyhKunfYencwyvpVtiN+6ylrsxfjA3bxO/mb6APl3bMn3iiXSMiXQ7JNPCWXI35ih9tHQjv52xgH7d2/HK1SfSPtom/zLus24ZY47C7CUbmPTqAvr3aM8rEy2xm8Bhyd2YIzRrcTG/nbGQIUkdeGXiCNq3scRuAod1yxhzBN5ZUMTtby5muKcTL1w5nJgo+yqZwGK/kcYcpje+W8dd7+RyUq/OPDch1dY4NQHJfiuNOQyvflPI3TOX8NM+XfjHFam0bhXudkjG+GR97sY00sv/zefumUs447h4S+wm4FnL3ZhGeP7LH3j4g+WcfXxXnho/lKgIS+wmsFnL3RgfsrOz8Xg8hIWF0blbAnc9+jTn9u/G0+OHWWI3QcGSuzEHyM7OJj09nYKCAlSVrZuKKfv4aX5CHpER9pUxwcF+U405QEZGxn5zrgNU7d3Nfffe41JExhw+S+7GHKC+Ze9sOTwTTCy5G3OAxKQkn+W2HJ4JJpbcjTnA0It/g0TsP12vLYdngk2jkruI5IvIEhFZJCI5TlknEZkjIqud545OuYjIkyKyRkRyRWRYU34AY/zpw9wNLI4ayNhbHrbl8ExQO5xx7qer6uY6+5OBuao6RUQmO/t3AecBfZzHicCzzrMxAW3d1gomv53LkKQOZF93G60eu8PtkIw5YkfTLTMGmOZsTwMuqlP+snp9DXQQke5H8XOMaXKV1TVMmrEQBP566VBahVuPpQlujf0NVuBjEZkvIulOWVdV3QDgPMc75QnAujqvLXLK9iMi6SKSIyI5paWlRxa9MX7yp3+tZPG67Tz6y0EkdYp2Oxxjjlpju2VOVtViEYkH5ojIigbqio8yPahANQvIAkhNTT3ouDHNZd7KEv7+xVrSTkxm1EA7yTShoVEtd1Utdp5LgJnACGBTbXeL81ziVC8C6o4lSwSK/RWwMf60acdubntjMcd1i+Xe0ce7HY4xfnPI5C4iMSISW7sNnAMsBWYBE5xqE4D3nO1ZwBXOqJmRQFlt940xgaS6Rrnl9UVU7K3mqfFDbZZHE1Ia0y3TFZgpIrX1X1XVj0TkO+ANEZkIFAKXOPVnA6OANUAFcJXfozbGD575bA1ffb+Fx8YOond8rNvhGONXh0zuqroWGOyjfAtwpo9yBW7wS3TGNJFvf9jK1E9WMWZIDy45IdHtcIzxOxvvZVqcbeV7uem1hSR3iibzFwNxzkqNCSm2WIdpUVSVO97KZfOPe3jnNyfT1ha2NiHKWu6mRXnpq3w+ydvE787rx8DE9m6HY0yTseRuWoyl68v4w+wVnNUvnqtO9rgdjjFNypK7aRF+3FPFpFcX0LltJH8cO9j62U3Isw5HE/JUlXtmLqFwawWvpZ9Ex5hIt0MypslZy92EvLfmF/HuomJuPqsvI3p2cjscY5qFJXcT0taU7OS+95ZxUq/O3HB6b7fDMabZWHI3IWt3ZTWTXl1IdGQ4j48bQniY9bOblsP63E3I+v2Hy1mxcScvXjWcru1aux2OMc3KWu4mJP1zyQamf11I+qm9OP3Y+EO/wJgQY8ndhJx1Wyu48+1cBid14PZzjnU7HGNcYcndhJTK6hpufG0hKDx16VAiI+xX3LRM1uduQsqfP17FwsLtPD1+mC2XZ1o0a9aYkPH5qlL+9vn3jD8xmfMH2XJ5pmWz5G5CQsmO3dz6+iKO7RrLfbZcnjHWLWOCX02NcusbiynfW8Vr40facnnGYC13EwKyvy3kyzWbuf+C/vTpasvlGQOW3E2QW799F1Nm5/HTPl0YNzzJ7XCMCRiNTu4iEi4iC0XkA2e/p4h8IyKrReR1EYl0yqOc/TXOcU/ThG5aOlUlY+YSFHjElsszZj+H03K/Ccirs/8oMFVV+wDbgIlO+URgm6r2BqY69Yzxu5kL1zNvZSl3/PxYG/ZozAEaldxFJBE4H3jO2RfgDOAtp8o04CJne4yzj3P8TLEmlfGz0p17eOiD5ZyQ0pErTvK4HY4xAaexLffHgTuBGme/M7BdVauc/SIgwdlOANYBOMfLnPr7EZF0EckRkZzS0tIjDN+0VA/MWkbFnmoe/eVAm+3RGB8OmdxFZDRQoqrz6xb7qKqNOPa/AtUsVU1V1dS4uLhGBWsMwEdLN/Lhkg3cdFYfesfHkp2djcfjISwsDI/HQ3Z2ttshGuO6xoxzPxm4UERGAa2Bdnhb8h1EJMJpnScCxU79IiAJKBKRCKA9sNXvkZsWqayiknvfW8rx3duRfmovsrOzSU9Pp6KiAoCCggLS09MBSEtLczNUY1x1yJa7qv5OVRNV1QOMAz5V1TTgM2CsU20C8J6zPcvZxzn+qaoe1HI35kj8/sPlbC3fy2NjB9EqPIyMjIx9ib1WRUUFGRkZLkVoTGA4mnHudwG3isgavH3qzzvlzwOdnfJbgclHF6IxXl+sKuXN+UVce2ovBiS0B6CwsNBn3frKjWkpDmv6AVWdB8xzttcCI3zU2Q1c4ofYjNmnfE8Vv3tnCb3iYrjxzD77ypOTkykoKDiofnJycnOGZ0zAsTtUTVD4479WUly2i8d+OWi/uWMyMzOJjt5/jHt0dDSZmZnNHaIxAcWSuwl4OflbmfbffCac5CHV02m/Y2lpaWRlZZGSkoKIkJKSQlZWll1MNS2eBMK1ztTUVM3JyXE7DBOAdldWM+rJf7OnsoaPbzmVmCibyNSYWiIyX1VTfR2zb4oJaE/OXc3a0nJevnqEJXZjDoN1y5iAtXR9GX//Yi1jT0jk1L52o5sxh8OSuwlIldU13PlWLp1iIrn3fFtZyZjDZee5JiBlfbGW5Rt28LfLTqB9dCu3wzEm6FjL3QScNSU/8sQnqxk1sBvnDujmdjjGBCVL7iagVNcod72dS5vIcB68cIDb4RgTtCy5m4Dy8n/zmV+wjftGH09cbJTb4RgTtCy5m4CxbmsFj320kp/1jePiYQmHfoExpl6W3E2zaWjedVXld+8sIUzgkYttPVRjjpaNljHN4lDzrr+ZU8SXazbz8Jj+JHRo42aoxoQEa7mbZtHQvOubduzm4Q+XM6JnJ9JOTHEpQmNCiyV30ywamnf93neXsreqhikXDyTM1kM1xi8suZtmUd/86l269eDj5Zu45ey+9Ipr28xRGRO6LLmbZuFr3vU2baKJOfkyBia059en9HQpMmNCk11QNc2idn71jIwMCgsLSU5Opt8F6ayOHcKjvxxERLi1M4zxJ/tGmWaTlpZGfn4+NTU1vPjRt+TFDOb6047h+B7t3A7NmJBzyOQuIq1F5FsRWSwiy0TkQae8p4h8IyKrReR1EYl0yqOc/TXOcU/TfgQTbHburuTumUvoE9+WG87o7XY4xoSkxrTc9wBnqOpgYAhwroiMBB4FpqpqH2AbMNGpPxHYpqq9galOPWP2mfLPFWzcsZtHxw4iKiL80C8wxhy2QyZ39frR2W3lPBQ4A3jLKZ8GXORsj3H2cY6fKXa7oXF8vXYL2d8UcvXJPRmW3NHtcIwJWY3qcxeRcBFZBJQAc4Dvge2qWuVUKQJqJwNJANYBOMfLgM4+3jNdRHJEJKe0tPToPoUJCjt3V3LHW4tJ6tSG287p63Y4xoS0RiV3Va1W1SFAIjAC6OermvPsq5V+0CrcqpqlqqmqmhoXZ0uotQT3vruU9dt2MfVXQ4iOtIFaxjSlwxoto6rbgXnASKCDiNR+QxOBYme7CEgCcI63B7b6I1gTvGYuLOLdRcXceGYfUj2d3A7HmJDXmNEycSLSwdluA5wF5AGfAWOdahOA95ztWc4+zvFPVfWglrtpOQq2lHPPzKUM93Rk0uk2OsaY5tCYc+PuwDQRCcf7x+ANVf1ARJYDr4nI74GFwPNO/eeBV0RkDd4W+7gmiNsEicrqGm58bRFhYcLj44bazUrGNJNDJndVzQWG+ihfi7f//cDy3cAlfonOBL2pc1axeN12nh4/zKbyNaYZWTPKNJmvvt/Ms59/z/+lJnH+oO5uh2NMi2LJ3TSJbeV7ufX1xfTsHMP9Fx7vdjjGtDiW3I3fqSp3vZ3LlvI9PHnpUBv2aIwLLLkbv8v+ppCPl2/irnOPY0BCe7fDMaZFsuRu/GrVpp08/MFyftqnC1efbHO0G+MWS+7Gb3ZXVnPjjIW0jYrgz78abEvmGeMi6ww1fjPlnytYsXEnL145nPjY1m6HY0yLZi134xefrtjES1/lc9XJHk4/Lt7tcIxp8Sy5m6NWsmM3t7+ZS7/u7Zh83nFuh2OMwZK7OUo1Ncptby6mYm8Vf710iC2+YUyAsORujspzX67l36s3c9/o/vSOj3U7HGOMw5K7OWJLisr4479Wcm7/blw6IsntcIwxdVhyN0ekfE8VN762kC5to5jyy4HYSorGBBYbCmmOyAOzlpG/pZwZ14ykQ3Sk2+EYYw5gLXdz2N5fXMyb84u44bTejOx10PK4xpgAYMndHJZ1Wyu4e+YShiZ34Kaz+rgdjjGmHpbcTaNVVddwy+uLUIUn/m8orWxVJWMClvW5m0b766dryCnYxhPjhpDcOdrtcIwxDWjMAtlJIvKZiOSJyDIRuckp7yQic0RktfPc0SkXEXlSRNaISK6IDGvqD2Ga3nf5W/nrp6u5eGgCY4YkuB2OMeYQGnNeXQXcpqr9gJHADSJyPDAZmKuqfYC5zj7AeUAf55EOPOv3qE2zKquo5ObXFpHUKZqHLhrgdjjGmEY4ZHJX1Q2qusDZ3gnkAQnAGGCaU20acJGzPQZ4Wb2+BjqIiC2gGaRUlbtnLmHTjt08MW4obaOsJ8+YYHBYV8RExAMMBb4BuqrqBvD+AQBqpwJMANbVeVmRU2aCRHZ2Nh6Ph7CwMOJ7JPH6a69y6zl9GZLUwe3QjDGN1OjkLiJtgbeBm1V1R0NVfZSpj/dLF5EcEckpLS1tbBimiWVnZ5Oenk5BQQGqyuaN69n+8dO0K/rG7dCMMYehUcldRFrhTezZqvqOU7yptrvFeS5xyouAuhONJALFB76nqmapaqqqpsbFxR1p/MbPMjIyqKio2K+seu9u7rknw6WIjDFHojGjZQR4HshT1b/UOTQLmOBsTwDeq1N+hTNqZiRQVtt9YwJfYWHhYZUbYwJTY1ruJwOXA2eIyCLnMQqYApwtIquBs519gNnAWmAN8A/gev+HbZpKUpLv2R2Tk5ObORJjzNE45NAHVf0S3/3oAGf6qK/ADUcZl3GBqjL04t+w7qkH0Ko9+8qjo6PJzMx0MTJjzOGy+8cN4E3sD76/nEVRAxl7y8MkJycjIqSkpJCVlUVaWprbIRpjDoMNWjaoKr//MI+Xvsrn16f0JOP8Uchjd7gdljHmKFjLvYVTVaZ8tILnv/yBK3/iIeP8frbwhjEhwJJ7C6aq/PnjVfz987VcNjKZ+y843hK7MSHCknsL9sTc1Tz12RouHZHEQxcOsMRuTAix5N5CPfXpah7/ZDWXnJBI5kUDCQuzxG5MKLHk3gL97fPv+dPHq7h4aAJTfjnIErsxIciSewvz3L/XMuWfK7hgcA/+eMlgwi2xGxOSLLm3IC/+5wd+/2Ee5w/sztRfWWI3JpRZcm8hXvm6gAffX87P+3fl8XFDiLD1T40JafYNbwFmfFvIve8u5ax+8fz10mG2sLUxLYB9y0PcGznruHvmEk4/No6n04YRGWH/5ca0BPZND2HvLCjirrdzOaV3F5697ASiIsLdDskY00wsuYeo9xat5/Y3F3NSr87844pUWreyxG5MS2LJPQR9mLuBW99YzHBPJ56bYIndmJbIknuI+WjpRm56bSHDkjvwwpXDiY60iT+NaYksuYeQT5Zv4rczFjAwsT0vXjWCmChL7Ma0VJbcQ8nP/y4AAAn9SURBVMRnK0u4PnsBx3dvx7SrR9DWErsxLZol9xDwxapSrn1lPn27teXliSfSrnUrt0MyxrjskMldRF4QkRIRWVqnrJOIzBGR1c5zR6dcRORJEVkjIrkiMqwpg2/pamqUl/+bzzUv53BMXFumTzyR9m0ssRtjGtdyfwk494CyycBcVe0DzHX2Ac4D+jiPdOBZ/4RpDrRq007G/u0rbn3kKdY/ezX/uuVnDDm+L9nZ2W6HZowJAIfsmFXVL0TEc0DxGOA0Z3saMA+4yyl/WVUV+FpEOohId1Xd4K+AW7o9VdU889n3PDNvDXtWfMGPnzzDnt27ACgoKCA9PR3AFrQ2poU70j73rrUJ23mOd8oTgHV16hU5ZcYP5hdsZfSTX/LE3NWMGtgd/W7GvsReq6KigoyMDJciNMYECn8PqfA1h6z6rCiSjrfrhuTkZD+HEVp27q7kj/9ayStfF9CjfRtevGo4px8bz1/Hr/NZv7CwsJkjNMYEmiNN7ptqu1tEpDtQ4pQXAUl16iUCxb7eQFWzgCyA1NRUn38ADMzN28Q97y5l447dTDjJwx0/P3bf+PXk5GQKCgoOeo39sTTGHGm3zCxggrM9AXivTvkVzqiZkUCZ9bcfmdKde5j06gImTsshtnUEb//mJzxwYf/9bkzKzMwkOjp6v9dFR0eTmZnZ3OEaYwLMIVvuIjID78XTLiJSBNwPTAHeEJGJQCFwiVN9NjAKWANUAFc1QcwhTVV5a34Rv/8wj117q7nt7L5c+7NjfE7VW3vRNCMjg8LCQpKTk8nMzLSLqcYYxDuwxV2pqamak5PjdhiuK9xSwe9m5vKfNVsY7unIHy4eRO/4tm6HZYwJUCIyX1VTfR2ze9QDQFV1DS/85wf+MmcVEWFh/P6iAYwfkUyYrXFqjDlCNv1AM8nOzsbj8RAWFobH49l3s9HS9WVc9Mx/eGT2Ck7pHcecW0/lspEpltiNMUfFWu7NIDs7m/T0dCoqKoD/3Ww0a3ExORH96RgdyTNpwzhvQDdELKkbY46e9bk3A4/H43PIYni7OG57/hPuHtWP9tE2J4wx5vBYn7vL6rupqGbnZh4dO6iZozHGtATW597EVJWuPXzPwGA3Gxljmoq13JtIwZZyZi5cz7sL11M59P+QkqfQyj37jtvNRsaYpmTJ3Y/KKir5YEkxMxesJ6dgGyJwUq/O3HD/TewcO5iH7r/XbjYyxjQLu6B6lPZW1fDZyhJmLljPpytK2FtdQ5/4tlw8LJExQ3rQo0Mbt0M0xoQou6DqZ6rKwnXbmblgPe/nFrO9opIubSO5bGQKFw9LoH+Pdjak0RjjKkvuB8jOzq53rpbCLRXefvRF6/lhczlREWGc078bFw9L4Ke9uxARbtenjTGBwZJ7Hb5uNrrmmnS+WrOFkq6pfJe/DfD2o//mtGM4b0A3Ym0xamNMALI+9zoautnotPvf5BdDE7hoaAIJ1o9ujAkA1ufegN2V1eRt2MHS4h0UNHCz0ZxbTrV+dGNM0GhRyf3HPVUsL97B0vVlLC0uY9n6Hawp/ZHqGu/ZS6t2cVSWlRz0uuTkZEvsxpigErTJvaELnwDbK/aybF8i38Gy9WX8sKWc2l6o+NgoBiS05+f9u9I/oT0DEtozb+Bf9utzB7vZyBgTnIIyufu68Dnx19fw8bKNtOl3GkuLyyjatmtf/YQObRiQ0I5fDE1gQEJ7+vdoR3y71ge9r61sZIwJFUF5QbWhC58/vfcN+vdox4CE9gzo4U3kHWMi/RmuMcYEhJC7oNrQLIuf3X5a8wZjjDEBqEnuuhGRc0VkpYisEZHJ/n7/+mZTtFkWjTHGy+/JXUTCgaeB84DjgUtF5Hh//ozMzEyio6P3K7MLn8YY8z9N0XIfAaxR1bWquhd4DRjjzx+QlpZGVlYWKSkpiAgpKSlkZWXZhU9jjHE0RZ97ArCuzn4RcOKBlUQkHUiHI+tOSUtLs2RujDH1aIqWu6+7fQ4akqOqWaqaqqqpcXFxTRCGMca0XE2R3IuApDr7iUBxE/wcY4wx9WiK5P4d0EdEeopIJDAOmNUEP8cYY0w9/N7nrqpVIjIJ+BcQDrygqsv8/XOMMcbUr0luYlLV2cDspnhvY4wxhxYQ0w+ISClw8HwC7uoCbHY7iMMQTPFarE0nmOINplghMONNUVWfI1ICIrkHIhHJqW/OhkAUTPFarE0nmOINplgh+OK1RT+NMSYEWXI3xpgQZMm9flluB3CYgilei7XpBFO8wRQrBFm81udujDEhyFruxhgTgiy5G2NMCLLkfgARSRKRz0QkT0SWichNbsd0KCISLiILReQDt2M5FBHpICJvicgK59/4JLdjqo+I3OL8DiwVkRkicvDCuy4SkRdEpEREltYp6yQic0RktfPc0c0Ya9UT6x+d34NcEZkpIh3cjLEuX/HWOXa7iKiIdHEjtsay5H6wKuA2Ve0HjARu8PdiI03gJiDP7SAa6QngI1U9DhhMgMYtIgnAjUCqqg7AO5XGOHejOshLwLkHlE0G5qpqH2Cusx8IXuLgWOcAA1R1ELAK+F1zB9WAlzg4XkQkCTgb8L3WZwCx5H4AVd2gqguc7Z14k0+Cu1HVT0QSgfOB59yO5VBEpB1wKvA8gKruVdXt7kbVoAigjYhEANEE2OymqvoFsPWA4jHANGd7GnBRswZVD1+xqurHqlrl7H6NdwbZgFDPvy3AVOBOfExjHmgsuTdARDzAUOAbdyNp0ON4f9lq3A6kEXoBpcCLTjfScyIS43ZQvqjqeuBPeFtoG4AyVf3Y3agapauqbgBvQwWIdzmexroa+KfbQTRERC4E1qvqYrdjaQxL7vUQkbbA28DNqrrD7Xh8EZHRQImqznc7lkaKAIYBz6rqUKCcwOk22I/TVz0G6An0AGJE5DJ3owpNIpKBtzs02+1Y6iMi0UAGcJ/bsTSWJXcfRKQV3sSerarvuB1PA04GLhSRfLxr1Z4hItPdDalBRUCRqtaeCb2FN9kHorOAH1S1VFUrgXeAn7gcU2NsEpHuAM5zicvxNEhEJgCjgTQN7JtujsH7h36x831LBBaISDdXo2qAJfcDiIjg7RPOU9W/uB1PQ1T1d6qaqKoevBf7PlXVgG1dqupGYJ2IHOsUnQksdzGkhhQCI0Uk2vmdOJMAvfh7gFnABGd7AvCei7E0SETOBe4CLlTVCrfjaYiqLlHVeFX1ON+3ImCY8zsdkCy5H+xk4HK8reBFzmOU20GFkN8C2SKSCwwBHnE5Hp+cs4u3gAXAErzflYC6/VxEZgD/BY4VkSIRmQhMAc4WkdV4R3VMcTPGWvXE+hQQC8xxvmd/czXIOuqJN6jY9APGGBOCrOVujDEhyJK7McaEIEvuxhgTgiy5G2NMCLLkbowxIciSuzHGhCBL7sYYE4L+HzvTkOErg7uaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "p0 = [700,85,0.7]\n",
    "def fun1(t,b1,b2,b3):\n",
    "    return b1 / (1 + b2*np.exp(-1*b3*t))\n",
    "y1 = np.empty (15,  dtype=np.float)\n",
    "b, _ = curve_fit(fun1, t, y, p0)\n",
    "print(\"[β1 β2 β3] =\", b)\n",
    "for i in range(0, 15):\n",
    "    y1[i] = fun1(t[i], b[0], b[1], b[2])\n",
    "plt.plot(t, y1, label='MNK')\n",
    "plt.plot(t, y, \"ko\",  label='Data')\n",
    "plt.legend()\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Построим модель вида  $y = \\frac{\\beta_1}{(1 + \\beta_2 e^{-\\beta_3 t})^{1/\\beta_4}}$ , применив МНК."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[β1 β2 β3] = [699.64135945 195.79500488   0.75962554   1.27922586]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deVxVdf7H8dcHXBBccEFTEbDcFwTFpX1xqRxTpyxTKq1+MU2WNdWUDmVlUbbM2PRrpl9UpiWpjWU65TRq5pQzWeIS7oEFihu4oXjZ+f7+uEcG5LIpcO69fJ6PB4977/d8z71vFD587/ee8z1ijEEppZR38bE7gFJKqdqnxV0ppbyQFnellPJCWtyVUsoLaXFXSikv1MjuAADt2rUzYWFhdsdQSimPsmnTpqPGmCBX29yiuIeFhZGYmGh3DKWU8igiklbRNp2WUUopL6TFXSmlvJAWd6WU8kJuMefuSkFBAenp6eTm5todRdWQn58fwcHBNG7c2O4oSjVYblvc09PTadGiBWFhYYiI3XFUNRljOHbsGOnp6XTt2tXuOEo1WG47LZObm0vbtm21sHsYEaFt27b6jstLJCQkEBYWho+PD2FhYSQkJNgdSVWT247cAS3sHkr/37xDQkICMTExOBwOANLS0oiJiQEgOjrazmiqGtx25K6UsldsbGxJYT/L4XAQGxtrUyLvUtfvirS4V0JEuPPOO0seFxYWEhQUxJgxYwCYP38+Pj4+JCUllfTp168fqampgPPkrKNHjwKwadMmunbtypYtW8q8xrp162jVqhWRkZH07t2b55577oIyP/vss7z22msAzJo1izVr1lTYd+vWraxcubLk8YoVK5gzZ84Fvb7yHvv27atRu6q+s++K0tLSMMaUvCuqzQKvxb0SAQEBbN++nZycHABWr15N586dy/QJDg4mLi6u0udJSkpiwoQJLFmyhMjIyHLbr7zySrZs2UJiYiILFy5k06ZNZbYXFhaeV/7Zs2czYsSICrefW9zHjh3LjBkzzuu1lPcJCQmpUbuqvvp4V6TFvQo33ngjX3zxBQCLFi1i0qRJZbaPGTOGHTt2sGfPHpf779q1i/Hjx/Phhx8yZMiQSl8rICCAQYMGsXfvXubPn8+tt97KTTfdxKhRowB49dVXGTx4MOHh4TzzzDMl+8XFxdGzZ09GjBhRJsfUqVNZunQpABs3buSyyy5jwIABDBkyhKysLGbNmsWSJUuIiIhgyZIlzJ8/nwcffBBwzq8OHz6c8PBwhg8fXjJamzp1KtOnT+eyyy7j4osvLnl+5X3i4uLw9/cv0+bv71/lYEZVrT7eFbn1B6pnPff3Hew8eKpWn7NPp5Y8c1PfKvvdfvvtzJ49mzFjxpCUlMQ999zDt99+W7Ldx8eHJ554ghdffJEFCxaU23/cuHEsXLiQK664osrXOnbsGBs2bODpp59m48aNfPfddyQlJdGmTRtWrVpFcnIyP/zwA8YYxo4dyzfffENAQACLFy9my5YtFBYWMnDgQAYNGlTmefPz85k4cSJLlixh8ODBnDp1Cn9/f2bPnk1iYiJvvvkm4JxmOuvBBx/krrvuYsqUKcybN4/p06fz2WefAXDo0CHWr1/P7t27GTt2LBMmTKjye1Oe5+yHprGxsezbt4+QkBDi4uL0w9RaEBISQlpa+WVhavNdUZUjdxHpKSJbS32dEpFHRKSNiKwWkWTrtrXVX0TkDRFJEZEkERlYa2ltEB4eTmpqKosWLWL06NEu+0yePJkNGzbwyy+/lNs2YsQI3n33XYqKiip8jW+//ZbIyEhGjRrFjBkz6NvX+Udn5MiRtGnTBoBVq1axatUqIiMjGThwILt37yY5OZlvv/2WX//61/j7+9OyZUvGjh1b7vn37NlDx44dGTx4MAAtW7akUaPK/65/9913TJ48GYA777yT9evXl2wbP348Pj4+9OnThyNHjlT6PMqzRUdHk5qaSnFxMampqbVS2PXwyvp5V1TlyN0YsweIABARX+AAsAyYAXxljJkjIjOsx08CNwLdra+hwFvW7Xmrzgi7Lo0dO5bHH3+cdevWcezYsXLbGzVqxGOPPcbLL79cbtubb77J/fffzwMPPMDbb7/t8vmvvPJKPv/883LtAQEBJfeNMcycOZPf/OY3Zfq8/vrrVR56aIy54MMTS+/ftGnTMs+tVHXp4ZVO9fGuqKZz7sOBvcaYNGAccHYeYgEw3ro/DvjAOG0AAkWkY62ktck999zDrFmz6N+/f4V9pk6dypo1a8jMzCzT7uPjw6JFi9izZw+zZs067wzXX3898+bNIzs7G4ADBw6QkZHBVVddxbJly8jJyeH06dP8/e9/L7dvr169OHjwIBs3bgTg9OnTFBYW0qJFC06fPu3y9S677DIWL14MOH8hqzOtpFRVGuLhlcYYCouKyS0oIjuvkCxHAUez8xhx0y2s37KTLEderb0rKq2mc+63A4us+x2MMYes8IdEpL3V3hnYX2qfdKvtUOknEpEYIAbc/9P34OBgHn744Ur7NGnShOnTp7vs17RpU5YvX87VV19Nhw4dmDZtWo0zjBo1il27dnHppZcC0Lx5cxYuXMjAgQOZOHEiERERhIaGcuWVV7rMtmTJEh566CFycnJo1qwZa9as4dprr2XOnDlEREQwc+bMMvu88cYb3HPPPbz66qsEBQXx/vvv1zizUueqyw8SExISyoyEZz//AuNvnUhuQRG5+cXkFhaRk19EbkEROQVF5BY4C27pxzkFReSVPC4ip1Sf/MJiCouN86uomKJiQ0HJraGo2FBYbPUpsu4XOftXJu7X/YgeGnrB3/+5pLpvq0WkCXAQ6GuMOSIiJ40xgaW2nzDGtBaRL4CXjDHrrfavgCeMMZtcPzNERUWZcy/WsWvXLnr37l3z70i5Bf3/U66EhYW5/CAxNDS05PyQc+UVFnHSUcDxM/mcOJPPcYd1e6aAE458jp/JZ/PaFXy34CWKC/JK9pNGTWlzw4M073ttjTL6NfbBr7EvzRr74lfy5UMTXx8a+QqNfHxo5CP/ve8r+PoIjX188PUVGvsIvj4+NLbaG/mW7v/ffc4+z8DQ1nRr37xGGUu+R5FNxpgoV9tqMnK/EdhsjDn7CdoREelojdo7AhlWezrQpdR+wTj/KCilGri4uLgyc+4ATfyacU30dJ7/fGfZ4u3I58SZArLzKj7Po6VfI9oENGHzJ2+VKewApjAPs3ERL8Q+VKpY+9CssS9NSz9u4otfI1+aNfGlia8PPj7esXxGTYr7JP47JQOwApgCzLFul5dqf1BEFuP8IDXr7PSNUqrhKC42HMzKITkjm70Z2aRkZJN8qittb3iIvDXzKDp1FN+W7Whx1V2sK+rJxh/20TqgCW0CmtDavwkXBzWntX8T2gQ0drb7NymzPdC/MY19nR8b+jyZ4TLDiYxD3DGs9qc8PEG1iruI+AMjgdKHaswBPhaRe4F9wK1W+0pgNJACOIC7ay2tUsrtFBYVk3bcQYpVwEt/5RT89xDgtgFN6Na+OXfeEU33R39D16DmBDVvSpsAZ6H2a+x73hnq47hxT1Ot4m6McQBtz2k7hvPomXP7GqDmnxgqpdxabkERvxw9Q7JVuPdmZJOccZrUow7yi4pL+nVq5ccl7ZszaUgI3do3L/lqE9CkzrK5mu5p6GfTesQZqkqp+pflKGDNriOs2XWEXYdOse+4g7MHfvgIhLTxp1v75lzXqwPd2jene/vmXNK+Oc2b1n9Z0bNpy9PirpQqcSw7j1U7j/CP7Yf5T8pRCosNF7X0Y2BoIGMjOpcU8a7tAi5oGqUuREdHN+hifi4t7pVIT09n2rRp7Ny5k+LiYsaMGcOrr75Kkybl314ePHiQ6dOnV7mQ1ujRo/noo48IDAystJ8rzz77LM2bN+fxxx8v1/7OO+8QFBTEmTNn6N+/Py+88AJ9+vSp9Pnmz5/PqFGj6NSpU42zKO9x5FQuX24/zD+2H+KHX45TbJyj8nuv6MoN/S5iQHCg1xxB0pB4zaqQtb1ehTGGm2++mfHjx5OcnMxPP/1Edna2yzPpCgsL6dSpU7VWSFy5cuV5Ffaq/O53v2Pr1q0kJyczceJErrvuunJny55r/vz5HDyoR6k2RPuPO3jnm5+5+a//ZuiLX/HMih0czc5n2rXd+GL6Ffzr99cwc3RvIkNaa2H3UF5R3Oti4fu1a9fi5+fH3Xc7D/bx9fVl7ty5zJs3D4fDUW5J3tTUVPr16wc4T6e+7bbbCA8PZ+LEiQwdOpSzJ2mdvYBHamoqvXv35r777qNv376MGjWqZN34d955h8GDBzNgwABuueWWcqdrV2XixImMGjWKjz76CHCu6z548GD69etHTEwMxhiWLl1KYmIi0dHRREREkJOT47Kf8h4/Z2bzl69TuOl/13PlK18Tt3IXuQXFPDayB2sevYo1j17NY6N60rdTK71UohfwiuJeF+tV7Nixo9zSuS1btiQkJISUlBTAuXLiggULWLt2bZl+f/3rX2ndujVJSUk8/fTT5S6+cVZycjLTpk1jx44dBAYG8sknnwBw8803s3HjRn788Ud69+7Ne++9V+P8Z1eOBOfyvRs3biy58Mjnn3/OhAkTiIqKIiEhga1bt9KsWTOX/ZRncPXO1RjD7sOnmLv6J66f+w3X/fFfvPrPPfj6CDNv7MW/fn8NKx++koeGd6db+xZ2fwuqlnnFnHtdrFdR0UqKpdtLL8lb2vr160vWmOnXrx/h4eEuX6Nr165EREQAMGjQoJLTr7dv385TTz3FyZMnyc7O5vrrrz+v/Gd9/fXXvPLKKzgcDo4fP07fvn256aabyu1T3X7KvbhaafHue/+HuC924Qi5FBEYHNaGZ27qw/V9L6JTYDObE6v64BXFvS5OYOjbt2/JSPqsU6dOsX//fi655BI2bdpUZkne0qo7nVF66VxfX9+SaZmpU6fy2WefMWDAAObPn8+6detqnH/Lli1ERUWRm5vLAw88QGJiIl26dOHZZ58lNze3XP/q9lPux9U714K8XNL++R7xK+9jVJ+LCGrRtIK9lbfyimmZulj4fvjw4TgcDj744AMAioqKeOyxx5g6dWq51zrXFVdcwccffwzAzp072bZtW41e+/Tp03Ts2JGCgoLz+tzgk08+YdWqVUyaNKmkQLdr147s7OwyH/qWXvK3sn7KvVX0DjXnxBGih4ZqYW+gvKK4R0dHEx8fT2hoKCJCaGgo8fHxF3TMq4iwbNky/va3v9G9e3d69OiBn58fL774YpX7PvDAA2RmZhIeHs7LL79MeHg4rVq1qvZrP//88wwdOpSRI0fSq1evau0zd+5cIiIi6N69OwsXLmTt2rUEBQURGBjIfffdR//+/Rk/fnzJ1ZjA+Q7h/vvvJyIigqZNm1bYT7mn9BMO7vsgEZ8W7Vxub8in3tcmj71ylDHG9q9BgwaZc+3cubNcm6coLCw0OTk5xhhjUlJSTGhoqMnLy7M5Vf3y5P8/d5dXUGT+8nWy6fnUStPrqX+Yu2P/ZPz9/Q1Q8uXv728WLlxod1SPt3DhQrf+twUSTQV11Svm3N2Nw+Hg2muvpaCgAGMMb731lssTn5Sqqf/sPcrTn21nb+YZru/bgVk39aVz4A0M791eT72vA5Udiefu/75a3OtAixYtOPfiI0pdiIzTubz4xS4+23qQLm2aMW9qFNf16lCyXU+9rxt1eeWouubWxd3UwoWdVf0zevJTrSkqNizckMZr/9xDXmEx04d354FrLnG7dV28lScvJey2H6j6+flx7NgxLRQexhjDsWPH8PPzszuKx9uy7wRj31zPMyt2EBESyD9/dxWPjuyhhb0e1cWRePXFbUfuwcHBpKenV7k+inI/fn5+BAcH2x3DY5105PPyl3tYvHEf7Vs05S+TBzK6/0X6LtYGnryUcLUvkF2XXF0gW6mGprjYsHRzOnP+sZusnALuviyMR0b2sGV9dOUZausC2UqpOrLr0Cme/mw7iWkniAptzfPj+9G7Y0u7YykPpsVdKRtl5xXy+uqfeP8/qbRq1phXJ4Rzy8BgXWZXXbDqXiA7EHgX6IfzQP57gD3AEiAMSAVuM8acEOfE4J9xXiTbAUw1xmyu9eRKeTBjDCu3HWb25zvIOJ3HpCEhPHF9TwL99XwIVTuqO3L/M/ClMWaCiDQB/IE/AF8ZY+aIyAxgBvAkcCPQ3foaCrxl3SqlcC4bMPPTbXybfJS+nVryf3cMIjKktd2xlJepsriLSEvgKmAqgDEmH8gXkXHANVa3BcA6nMV9HPCBdWrsBhEJFJGOxphDtZ5eKQ+z/7iDiW9/x+ncQp4b25c7hoXiq1Mwqg5U5zj3i4FM4H0R2SIi74pIANDhbMG2bttb/TsD+0vtn261lSEiMSKSKCKJerijagjSTzgY/tsXSHx5Mjuev5FnJl/D4kUf2R1LeanqFPdGwEDgLWNMJHAG5xRMRVwNQ8odb2mMiTfGRBljooKCgqoVVilPdfBkDiN+G0fKsj9RkJVRa5eDVKoi1Snu6UC6MeZ76/FSnMX+iIh0BLBuM0r171Jq/2BAr8KsGqzDWblMfmcDP3/5LqYgr8y2C70cpFIVqbK4G2MOA/tFpKfVNBzYCawAplhtU4Dl1v0VwF3iNAzI0vl21VBlnHIW9qPZ+RSecj396AmLUCnPU92jZR4CEqwjZX4G7sb5h+FjEbkX2AfcavVdifMwyBSch0LeXauJlfIQmafzmPTOBg6fyuWDe4YwYYHnLkKlPE+1irsxZivg6hTX4S76GmDaBeZSyqMdzc5j8jsbOHgyl/l3DyYqrA1xcXFlLmQNnrMIlfI8brsqpFKe6viZfO5493v2n3Awb+pghl7cFqiby0EqVRFdOEypWnTSkc/kd75nb2Y286YO5vJurq9vqlRt0IXDlKoHWY4C7njve1Iys3n3rigt7MpWOi2jVC04lVvAXfO+56fD2bx9xyCu6qHnbih7aXFX6gKdzi1gyrwf2HnoFH+NHsi1vdpXvZNSdUyLu1IXIDuvkKnvb2RbehZvTh7IiD4dqt5JqXqgc+5KnSdHfiH3vL+RrftP8uakSK7ve5HdkZQqoSN3pc5DTn4R985PJDHtOK9PjODG/h3tjqRUGVrclaqh3IIi7vsgke9/OcbciRHcNKCT3ZGUKkenZZSqgdyCImI+3MS/9x7ltQkDGBdRbjVrpdyCjtyVqqa8wiJ+u3AT3/yUycu3hHPLoGC7IylVIS3uSlVDfmEx0xI28/WeTF66uT+3RXWpeielbKTFXakqFBQV89CizazZlcEL4/sxaYiu4qjcnxZ3pVxISEggLCwMHx8f2nTozCcfLy655qlSnkCLu1LnSEhIICYmhrS0NIwxZB87zOnVf6XRL/+2O5pS1abFXalzxMbGlllzHSA/L0cvh6c8ihZ3pc5R0WXv9HJ4ypNocVfqHMHBro+E0cvhKU+ixV2pc/Qbdz/SqGmZNr0cnvI01SruIpIqIttEZKuIJFptbURktYgkW7etrXYRkTdEJEVEkkRkYF1+A0rVpi+SDrEzIJxJj72gl8NTHq0myw9ca4w5WurxDOArY8wcEZlhPX4SuBHobn0NBd6ybpVya4eycvjDsm0M6BLI/PsfpfGcx+2OpNR5u5BpmXHAAuv+AmB8qfYPjNMGIFBEdMk85daKiw2PLvmRgqJi/jwxgsa+OmOpPFt1f4INsEpENolIjNXWwRhzCMC6PXv5mc7A/lL7plttZYhIjIgkikhiZmbm+aVXqpa88+3PfPfzMZ69qS9h7QLsjqPUBavutMzlxpiDItIeWC0iuyvpKy7aTLkGY+KBeICoqKhy25WqL9sPZPHaqj3c0Pcibo3SxcCUd6jWyN0Yc9C6zQCWAUOAI2enW6zbDKt7OlD6WLJg4GBtBVaqNuXkF/HIkq20CWjCSzf3R8TV2EQpz1NlcReRABFpcfY+MArYDqwApljdpgDLrfsrgLuso2aGAVlnp2+Ucjcv/WMXKRnZvHbrAFoHNLE7jlK1pjrTMh2AZdaIphHwkTHmSxHZCHwsIvcC+4Bbrf4rgdFACuAA7q711ErVgrW7j/DBd2n8zxVdubJ7kN1xlKpVVRZ3Y8zPwAAX7ceA4S7aDTCtVtIpVUeOZufxxNIkel3Ugt/f0NPuOErVOr3MnmpwjDE8sTSJU7mFfHTfMJo28rU7klK1Tg/mVQ3Owg1prN2dwcwbe9GjQwu74yhVJ7S4qwYlJeM0L3yxi6t7BDH1sjC74yhVZ7S4qwYjv7CYhxdvJaBpI169NVwPe1ReTefcVYPxx9V72HHwFO/cFUX7Fn52x1GqTunIXTUI/9l7lPhvfmby0BBG9ulgdxyl6pwWd+X1shwFPPbxj3RtG8BTv+ptdxyl6oVOyyivZozhD59tI/N0HsseuBz/JvojrxoGHbkrr/bp5gN8kXSI343sQf/gVnbHUareaHFXXmvfMQezlm9nSNc23H/1JXbHUapeaXFXXqmwqJhHlmzBx0eYOzECXx897FE1LDoBqbzSX77ey+Z9J/nz7RF0Dmxmdxyl6p2O3JXX2bzvBG+sTWZ8RCfGRZS7CJhSDYIWd+VVsvMKeWTxVi5q6cfs8f3sjqOUbXRaRnmV51bsIP2EgyW/uZSWfo3tjqOUbXTkrrzGP7Yd4m+b0nngmm4MDmtjdxylbKXFXXmFQ1k5zPh0GwOCW/HwiO52x1HKdlrclccrLjY8/rcfyS8s5vXbI2nsqz/WSulvgfJ4izfu598px3h6TB+6tguwO45SbqHaxV1EfEVki4h8bj3uKiLfi0iyiCwRkSZWe1PrcYq1PaxuoisFh7NyeWnlLi67pC2ThnSxO45SbqMmI/eHgV2lHr8MzDXGdAdOAPda7fcCJ4wx3YC5Vj+lap0xhqeXb6eguJiXbu6vF99QqpRqFXcRCQZ+BbxrPRbgOmCp1WUBMN66P856jLV9uOhvnaoDK7cdZvXOIzw6sgehbXU6RqnSqjtyfx14Aii2HrcFThpjCq3H6cDZUwE7A/sBrO1ZVv8yRCRGRBJFJDEzM/M846uG6qQjn2dWbKd/51bcc3lXu+Mo5XaqLO4iMgbIMMZsKt3soqupxrb/NhgTb4yJMsZEBQUFVSusUme98MUuTjoKePmWcJYsXkRYWBg+Pj6EhYWRkJBgdzylbFedM1QvB8aKyGjAD2iJcyQfKCKNrNF5MHDQ6p8OdAHSRaQR0Ao4XuvJVYP1bXImSzelM+3aS9jy9d+JiYnB4XAAkJaWRkxMDADR0dF2xlTKVlWO3I0xM40xwcaYMOB2YK0xJhr4GphgdZsCLLfur7AeY21fa4wpN3JX6nw48guZ+ek2Lg4K4KHruhMbG1tS2Ev6OBzExsbalFAp93Ahx7k/CTwqIik459Tfs9rfA9pa7Y8CMy4solL/9cdVP5F+Ioc5N4fj19iXffv2uexXUbtSDUWNFg4zxqwD1ln3fwaGuOiTC9xaC9mUKmPr/pO8/+9fuGNYCEO6OteOCQkJIS0trVzfkJCQ+o6nlFvRM1SVR8gvLObJpUm0b+HHkzf0KmmPi4vD39+/TF9/f3/i4uLqO6JSbkWLu/II//evvew5cpoXxvejRamlfKOjo4mPjyc0NBQRITQ0lPj4eP0wVTV4up67cnspGad5c20KY8I7MqJPh3Lbo6OjtZgrdQ4duSu3VlxsePKTbfg39eXZsX3tjqOUx9DirtzahxvS2JR2gqd/1Yd2zZvaHUcpj6HFXbmtAydzeOXL3VzZvR03D9QLXStVE1rclVsyxhC7bBsGePHXuuKjUjWlxV25peVbD7JuTyaPj+pJlzb+Ve+glCpDi7tyO8ey83ju7zuI6BLIlMvC7I6jlEfS4q7czuzPd5KdV8grE8Lx9dHpGKXOhxZ35VbW7j7C8q0HeeCabvTo0MLuOEp5LC3uqt4kJCRUuu56dl4hTy3bTvf2zXng2ktsSqmUd9AzVFW9SEhIqHLd9Ve+3M2hU7ksvf8ymjbytS2rUt5AR+6qXlS17npi6nE+3JDGlEvDGBTa2o6ISnkVLe6qXlS27npuQRFPfpJEp1bN+P31Pes5mVLeSYu7qhcVra8eEhLCX75OYW/mGV68uT8BTXWmUKnaoMVd1YuK1l2f9sRTvLVuLzdHdubqHnqhdKVqixZ3VS9crbv+f2/Hs970plWzxjw9po/dEZXyKlrcVb2Jjo4mNTWV4uJiUlNTyQu9lB/3n+SZsX1pHdDE7nhKeZUqi7uI+InIDyLyo4jsEJHnrPauIvK9iCSLyBIRaWK1N7Uep1jbw+r2W1CeaP9xB39c9RPDe7XnpvCOdsdRyutUZ+SeB1xnjBkARAA3iMgw4GVgrjGmO3ACuNfqfy9wwhjTDZhr9VOqhDGGmZ9uw9dHeOHX/XTFR6XqQJXF3ThlWw8bW18GuA5YarUvAMZb98dZj7G2Dxf97VWlLN2UzvqUozx5Yy86tmpmdxylvFK15txFxFdEtgIZwGpgL3DSGFNodUkHzl5NoTOwH8DangW0dfGcMSKSKCKJmZmZF/ZdKI+ReTqPF77YxZCwNkQPcX14pFLqwlWruBtjiowxEUAwMATo7aqbdetqlG7KNRgTb4yJMsZEBQXpIXANgTGGpz7bRk5BES/d0h8fXfFRqTpTo6NljDEngXXAMCBQRM6ecRIMHLTupwNdAKztrYDjtRFWebaFG9L4544jPD6qB5cENbc7jlJerTpHywSJSKB1vxkwAtgFfA1MsLpNAZZb91dYj7G2rzXGlBu5q4Zlx8Esnv98F9f2DOJ/rrjY7jhKeb3qnOvdEVggIr44/xh8bIz5XER2AotF5AVgC/Ce1f894EMRScE5Yr+9DnIrD5KdV8iDH22hTUAT/nhbhE7HKFUPqizuxpgkINJF+88459/Pbc8Fbq2VdMrjGWP4w6fbSDt2hsUxl9JGT1ZSql7oGaqqTn2cuJ8VPx7k0ZE9GNK1jd1xlGowtLirOrPn8GmeWbGDK7q147fXdLM7jlINihZ3VScc+YVM+2gzzZs2Zu7ECL3QtVL1TBfPVnVi1vId7M3MZuG9Qwlq0dTuOEo1ODpyV7Xuk03pLN2UzlGbQIMAAA9jSURBVEPXduPybu3sjqNUg6TFXdWqlIxsnl6+naFd2/DwiB52x1GqwdLirmpNbkERD360Gb/Gvvz59kidZ1fKRjrnrmrN7M93svvwaebfPZiLWvnZHUepBk1H7qpW/P3Hg3z0/T7uv/oSrunZ3u44SjV4WtzVBUs9eoaZn25jUGhrHhul8+xKuQMt7uqC5BUW8eCizfj6CG9MiqSxr/5IKeUOdM5dXZCXVu5m+4FTvHNXFJ0D9apKSrkLHWap8/bl9sPM/08q91zelZF9OtgdRylVihZ3dV72H3fwxNIfCQ9uxYwbe9kdRyl1Di3uqsbyC4t5aNEWjIE3Jw2kSSP9MVLK3eicu6qx11btYev+k/w1eiAhbf3tjqOUckGHXKpG1u4+Qvw3P3PnsFBG9+9odxylVAW0uKtqO5SVw6Mf/0jvji2J/VVvu+MopSpRnQtkdxGRr0Vkl4jsEJGHrfY2IrJaRJKt29ZWu4jIGyKSIiJJIjKwrr8JVfcKi4qZvmgLBYXF/GVyJH6Nfe2OpJSqRHVG7oXAY8aY3sAwYJqI9AFmAF8ZY7oDX1mPAW4EultfMcBbtZ5a1bu5a35iY+oJXry5PxcHNbc7jlKqClUWd2PMIWPMZuv+aWAX0BkYByywui0Axlv3xwEfGKcNQKCI6OSsB/vmp0z+um4vE6O6MC6is91xlFLVUKM5dxEJAyKB74EOxphD4PwDAJxdLaozsL/UbulWm/IQCQkJhIWF4ePjQ5eQEKbE/pHu7Zvz7Ni+dkdTSlVTtYu7iDQHPgEeMcacqqyrizbj4vliRCRRRBIzMzOrG0PVsYSEBGJiYkhLS8MYQ/r+/aQtf53rm6bQrInOsyvlKapV3EWkMc7CnmCM+dRqPnJ2usW6zbDa04EupXYPBg6e+5zGmHhjTJQxJiooKOh886taFhsbi8PhKNNmCvL431eetymRUup8VOdoGQHeA3YZY/5UatMKYIp1fwqwvFT7XdZRM8OArLPTN8r97du3r0btSin3VJ2R++XAncB1IrLV+hoNzAFGikgyMNJ6DLAS+BlIAd4BHqj92KquhISE1KhdKeWeqlx+wBizHtfz6ADDXfQ3wLQLzKVscvtvn+CVpx7FFOaVtPn7+xMXF2djKqVUTekZqqrE50kHWZoVRmT0k3TpEoKIEBoaSnx8PNHR0XbHU0rVgC4cpgD4ZFM6v1/6I4NCWzNv1lO08HvO7khKqQugxV3x0ff7+MOybVzerS3v3BWFfxP9sVDK0+lvcQP3/r9/4bm/7+TankG8dccgXTNGKS+hxb0Be2vdXl7+cjfX9+3A/+pFN5TyKlrcGyBjDK+vSebPXyVz04BO/Om2ATT21cKulDfR4t7AGGOY8+Vu3v7Xz0wYFMzLt4Tj61PRka5KKU+lxb0BKS42zP58J/P/k8odw0KYPbYfPlrYlfJKWtwbiOJiQ+xn21j0w37uvaIrT/2qN86VJZRS3kiLewNQWFTME0uT+HTLAaZdewmPj+qphV0pL6fF3csVFBXzyOKtfLHtEI+N7MFDw7vbHUkpVQ+0uHuxvMIipiVsYc2uI8SO7s19V11sdySlVD3R4u6lcvKL+M3CTXzzUyazx/XlrkvD7I6klKpHWty90Jm8Qu5dsJHvfznOK7eEc9vgLlXvpJTyKlrcvcyp3AKmzvuBH9OzmHtbBOMj9fK1SjVEWty9yElHPnfN+4GdB0/x5qRIbuzf0e5ISimbaHH3Ekez87jj3e/5OfMMb985iOG9O9gdSSllIy3uXuDIqVwmv7OBAydzeG9qFFd21wuOK9XQaXH3cNvSs3hw0WaOns5j/t1DGHZxW7sjKaXcQJVLAYrIPBHJEJHtpdraiMhqEUm2bltb7SIib4hIiogkicjAugzfkJ3KLeCZ5dsZ95f1nMkr4sP/GaqFXSlVojrrvM4HbjinbQbwlTGmO/CV9RjgRqC79RUDvFU7MdVZxhiWbz3Ada/9iw83pBGRt52M+HuJCmtLWFgYCQkJdkdUSrmBKou7MeYb4Pg5zeOABdb9BcD4Uu0fGKcNQKCI6CEbtSQlI5vod7/n4cVb6RzoR0ynQ/zz7dns378PYwxpaWnExMRogVdKnfecewdjzCEAY8whEWlvtXcG9pfql261HTr/iConv4g3v04m/puf8Wvsy/Pj+zF5SAiXXHwHDoejTF+Hw0FsbCzR0dE2pVVKuYPa/kDV1VKDxmVHkRicUzeEhITUcgzvsXb3EWYt30H6iRxujuzMzNG9CWrRFIB9+/a53KeidqVUw3G+xf2IiHS0Ru0dgQyrPR0ofa57MHDQ1RMYY+KBeICoqCiXfwAasgMnc3huxQ5W7TxCt/bNWRwzrNwHpiEhIaSlpZXbV/9YKqXO98KZK4Ap1v0pwPJS7XdZR80MA7LOTt+o6ikoKubtf+1lxB//xTfJmTx5Qy9WTr/S5ZEwcXFx+Pv7l2nz9/cnLi6uvuIqpdxUlSN3EVkEXAO0E5F04BlgDvCxiNwL7ANutbqvBEYDKYADuLsOMnutH345zlOfbeOnI9mM6N2BZ8f2Ibi1f4X9z86rx8bGsm/fPkJCQoiLi9P5dqUUYoz9MyJRUVEmMTHR7hi2OZadx4srd/PJ5nQ6Bzbj2bF9GdlHlw9QSlVORDYZY6JcbdMzVG1UXGxYtHEfr3y5hzN5hfz2mkt46Lpu+DfR/xal1IXRKlJPEhISykyf/Pbxp/hOerN1/0mGXdyG58f1o3uHFnbHVEp5CS3u9SAhIYGYmJiSY9LT0tKY+bsHCR3/CHP/8CDjIzrrBauVUrXqfI+WUTUQGxtb7mQjU5hH0feL+HVksBZ2pVSt05F7HSouNmxMPU5aBScVpafvd9mulFIXSot7HUjJyGbZlnQ+23KQAydzaNQyiMKsjHL99GQjpVRd0eJeS45m57Fi60E+23qApPQsfASu6B7E76/vyYner/HgA/eXmZrRk42UUnVJi/sFyMkvYvWuIyzbnM43yUcpKjb07dSSp37Vm7EDOtG+pZ+zY+SdNGnkoycbKaXqjZ7EVENFxYbvfz7Gp1sO8OX2w2TnFdKplR/jIjvz68jO9NDDGZVS9URPYqoFew6f5tMt6SzfcpDDp3Jp3rQRo/tfxPjIzgzr2hYfHz3iRSnlPrS4n6P0yUadg4O5fsrvONAuil2HTtHIR7i6RxCxv+rNyD4d8Gvsa3dcpZRySadlSlnw4Yfc/5vfkJuTU9ImjZoSGf0k02PuZkx4R9o2b2pjQqWU+i+dlnGhqNiQkpHNj+knSUo/SVJ6Fv+IfZzCUoUdnCcbHVu3gCnzn7MpqVJK1VyDKO7GGNKOOaxCnkVS+km2HzhFTkERAC2aNqJf51YUns50ub9e2Ugp5Wm8srgfzsotMyJPSs8iK6cAgKaNfOjbqSUTB3chPLgV4cGBXNwuAB8f4bsX9cpGSinv4LHFvfQHn+07dmbUlEdo0vMqktKzyDidB4Cvj9CzQwtG97+I8OBAwoNb0aNDCxr7ul5SJy4urswCX6AnGymlPJNHFvdzV1k8cjCdha/G0vvWxxk1bkLJiLxvp5Y1OqJFr2yklPIWHnm0TFhYmMvpk9DQUFJTU2sxmVJKua/KjpbxyCV/K/qAUz/4VEoppzop7iJyg4jsEZEUEZlR289f0Qec+sGnUko51XpxFxFf4C/AjUAfYJKI9KnN14iLi8Pf379Mm37wqZRS/1UXI/chQIox5mdjTD6wGBhXmy8QHR1NfHw8oaGhiAihoaHEx8frB59KKWWpi6NlOgOlLzGUDgw9t5OIxAAxcH7TKdHR0VrMlVKqAnUxcne1PGK5Q3KMMfHGmChjTFRQUFAdxFBKqYarLop7OtCl1ONg4GAdvI5SSqkK1EVx3wh0F5GuItIEuB1YUQevo5RSqgK1PudujCkUkQeBfwK+wDxjzI7afh2llFIVq5PlB4wxK4GVdfHcSimlquYWyw+ISCZQfj0Be7UDjtodogY8Ka9mrTuelNeTsoJ75g01xrg8IsUtirs7EpHEitZscEeelFez1h1PyutJWcHz8nrk2jJKKaUqp8VdKaW8kBb3isXbHaCGPCmvZq07npTXk7KCh+XVOXellPJCOnJXSikvpMVdKaW8kBb3c4hIFxH5WkR2icgOEXnY7kxVERFfEdkiIp/bnaUqIhIoIktFZLf1b3yp3ZkqIiK/s34GtovIIhHxsztTaSIyT0QyRGR7qbY2IrJaRJKt29Z2ZjyrgqyvWj8HSSKyTEQC7cxYmqu8pbY9LiJGRNrZka26tLiXVwg8ZozpDQwDptX2xUbqwMPALrtDVNOfgS+NMb2AAbhpbhHpDEwHoowx/XAupXG7vanKmQ/ccE7bDOArY0x34CvrsTuYT/msq4F+xphw4CdgZn2HqsR8yudFRLoAIwG3v6anFvdzGGMOGWM2W/dP4yw+ne1NVTERCQZ+Bbxrd5aqiEhL4CrgPQBjTL4x5qS9qSrVCGgmIo0Af9xsdVNjzDfA8XOaxwELrPsLgPH1GqoCrrIaY1YZYwqthxtwriDrFir4twWYCzyBi2XM3Y0W90qISBgQCXxvb5JKvY7zh63Y7iDVcDGQCbxvTSO9KyIBdodyxRhzAHgN5wjtEJBljFllb6pq6WCMOQTOgQrQ3uY81XUP8A+7Q1RGRMYCB4wxP9qdpTq0uFdARJoDnwCPGGNO2Z3HFREZA2QYYzbZnaWaGgEDgbeMMZHAGdxn2qAMa656HNAV6AQEiMgd9qbyTiISi3M6NMHuLBUREX8gFphld5bq0uLugog0xlnYE4wxn9qdpxKXA2NFJBXntWqvE5GF9kaqVDqQbow5+05oKc5i745GAL8YYzKNMQXAp8BlNmeqjiMi0hHAus2wOU+lRGQKMAaINu590s0lOP/Q/2j9vgUDm0XkIltTVUKL+zlERHDOCe8yxvzJ7jyVMcbMNMYEG2PCcH7Yt9YY47ajS2PMYWC/iPS0moYDO22MVJl9wDAR8bd+Jobjph/+nmMFMMW6PwVYbmOWSonIDcCTwFhjjMPuPJUxxmwzxrQ3xoRZv2/pwEDrZ9otaXEv73LgTpyj4K3W12i7Q3mRh4AEEUkCIoAXbc7jkvXuYimwGdiG83fFrU4/F5FFwHdATxFJF5F7gTnASBFJxnlUxxw7M55VQdY3gRbAauv37P9sDVlKBXk9ii4/oJRSXkhH7kop5YW0uCullBfS4q6UUl5Ii7tSSnkhLe5KKeWFtLgrpZQX0uKulFJe6P8BahXZKaSR9EkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "p0 = [700,200,0.7,1.3]\n",
    "def fun2(t,b1,b2,b3,b4):\n",
    "    return b1 / ((b2*np.exp(-1*b3*t) + 1)**(1/b4))\n",
    "y2 = np.empty (15,  dtype=np.float)\n",
    "b, _ = curve_fit(fun2, t, y, p0)\n",
    "print(\"[β1 β2 β3] =\", b)\n",
    "for i in range(0, 15):\n",
    "    #y2[i] = fun2(t[i], 700, 70, 0.01, 0.01)\n",
    "    y2[i] = fun2(t[i], b[0], b[1], b[2], b[3])\n",
    "plt.plot(t, y2,  label='MNK Prediction')\n",
    "plt.plot(t, y, \"ko\",  label='Original Data')\n",
    "plt.legend()\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Общую ошибку обеих моделей вычислим по данной формуле: $\\frac{\\sum_{i=1}^n (y_i - f_i(x))^2}{n} $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ошибка модели 1 =  595.3255315004633\n",
      "Ошибка модели 2 =  585.760327350927\n"
     ]
    }
   ],
   "source": [
    "er1 = np.array (y - y1)\n",
    "for i in range (0,15):\n",
    "    er1[i] = er1[i]**2\n",
    "error1 = sum(er1)/15\n",
    "er2 = np.array (y - y2)\n",
    "for i in range (0,15):\n",
    "    er2[i] = er2[i]**2\n",
    "error2 = sum(er2)/15\n",
    "print(\"Ошибка модели 1 = \", error1)\n",
    "print(\"Ошибка модели 2 = \", error2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Вывод:\n",
    "Вторая модель оказалась немного точнее первой. Очевидно, что точность решения напрямую зависит от количества переменных, по которым составляется уравнение."
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}