SimplifiedSIR.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simplified SIR model with vaccinations\n",
    "\n",
    "The standard SIR model reads:\n",
    "\n",
    "$\\frac{dS}{dt} = - \\beta I \\frac{S}{N}$\n",
    "\n",
    "$\\frac{dI}{dt} =   \\beta I \\frac{S}{N} - \\gamma I$\n",
    "\n",
    "$\\frac{dR}{dt} =  \\gamma I$\n",
    "\n",
    "I use the following approximation:\n",
    "\n",
    "$\\frac{S}{N} \\sim 1$\n",
    "\n",
    "Hence, I can use only the second and the third equation.\n",
    "In addition, one can introduce vaccinations too.\n",
    "Another approximation is that the vaccination rate is constant in time, which is not hard to believe, and mostly true, looking at the available data.\n",
    "With this, I write:\n",
    "\n",
    "$\\frac{dI}{dt} =   \\beta I \\frac{S}{N} - \\gamma I - \\alpha$\n",
    "\n",
    "$\\frac{dR}{dt} =  \\gamma I$\n",
    "\n",
    "where $\\alpha$ comes from data.\n",
    "There are two scenarios here:\n",
    "\n",
    "$\\alpha_{IT} \\sim 0.2$\n",
    "\n",
    "$\\alpha_{UK} \\sim 0.6$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.optimize as scopt\n",
    "import scipy.integrate as scint\n",
    "import scipy.optimize as scopt\n",
    "\n",
    "def pend(y, t, b, c):\n",
    "    theta, omega = y\n",
    "    dydt = [omega, -b*omega - c*np.sin(theta)]\n",
    "    return dydt\n",
    "\n",
    "\n",
    "def test_pend(y, t, b, c):\n",
    "    theta, omega = y\n",
    "    dydt = [omega, -b*omega - c*np.sin(theta)]\n",
    "    return dydt\n",
    "\n",
    "\n",
    "def SimplifiedSIR(y, t, beta, tau, gamma):\n",
    "    I, R = y\n",
    "    dydt = [ (beta*np.exp(-t/tau) - gamma)*I, gamma*I]\n",
    "    return dydt\n",
    "\n",
    "def SimplifiedSIR_2(y, t, beta, gamma):\n",
    "    I, R = y\n",
    "    dydt = [ (beta - gamma)*I, gamma*I]\n",
    "    return dydt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Taking the constants from the following paper:\n",
    "https://arxiv.org/pdf/2005.08724.pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "beta_0  = 0.384\n",
    "tau_0   = 26.33\n",
    "gamma_0 = 0.048\n",
    "\n",
    "y_0 = [1.,1.]\n",
    "\n",
    "\n",
    "guess = [0.5, 30, 0.1]\n",
    "t     = np.linspace(0, 1, 11)\n",
    "t2    = np.linspace(0, 1, 14)\n",
    "\n",
    "sol   = scint.odeint(SimplifiedSIR, y_0, t, args=(beta_0, tau_0, gamma_0))\n",
    "\n",
    "# popt, pcov = scopt.curve_fit(test_pend, guess, t, sol)\n",
    "# popt, pcov = scopt.curve_fit(SimplifiedSIR, guess, t, sol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.        , 1.        ],\n",
       "       [1.03409552, 1.00488143],\n",
       "       [1.06919821, 1.00992892],\n",
       "       [1.10533249, 1.01514738],\n",
       "       [1.14252319, 1.02054181],\n",
       "       [1.1807956 , 1.02611733],\n",
       "       [1.22017546, 1.03187922],\n",
       "       [1.26068891, 1.03783283],\n",
       "       [1.30236256, 1.04398369],\n",
       "       [1.34522346, 1.05033741],\n",
       "       [1.38929909, 1.05689978]])"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "# popt, pcov = scopt.curve_fit(SimplifiedSIR, guess, t, sol)\n",
    "\n",
    "\n",
    "def f(t, a0, a1, beta, tau, gamma): \n",
    "    y_00 = [a0, a1]\n",
    "    tspan = np.hstack([[0],np.hstack([t])])\n",
    "    return scint.odeint(SimplifiedSIR, y_00, tspan, args=(beta,tau,gamma))[1:,0]\n",
    "#     return scint.odeint(SimplifiedSIR, [a0,a1], tspan, args=(beta,tau,gamma))[1:,0]\n",
    "\n",
    "def f0(t, beta,tau,gamma): \n",
    "    tspan = np.hstack([[0],np.hstack([t])])\n",
    "    return scint.odeint(SimplifiedSIR, y_0, tspan, args=(beta,tau,gamma))[1:,0]\n",
    "\n",
    "\n",
    "guess2 = [1., 1., 0.5, 30, 0.1]\n",
    "\n",
    "popt, pcov = scopt.curve_fit(f, t, sol[:,0], p0=guess2)\n",
    "# popt, pcov = scopt.curve_fit(f0, t, sol[:,0], p0=guess)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.00000083,  1.33371254,  0.48405754, 33.32701271,  0.1480694 ])"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "popt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Loading the data\n",
    "\n",
    "Here I am writing by hand (sigh) the data I want to analyse..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[255.28,   4.71],\n",
       "       [264.7 ,   4.78],\n",
       "       [274.13,   4.68],\n",
       "       [282.37,   4.61],\n",
       "       [291.26,   4.58],\n",
       "       [301.75,   4.65],\n",
       "       [308.78,   4.73],\n",
       "       [317.1 ,   4.83],\n",
       "       [328.8 ,   4.89],\n",
       "       [336.03,   4.93],\n",
       "       [337.87,   5.1 ],\n",
       "       [344.16,   5.18],\n",
       "       [347.55,   5.14],\n",
       "       [354.17,   5.22]])"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_IT = np.array([[255.28,4.71], [264.70,4.78], [274.13,4.68], [282.37,4.61], [291.26,4.58], [301.75,4.65], [308.78,4.73],  #26 Feb\n",
    "                    [317.10,4.83], [328.80,4.89], [336.03,4.93], [337.87,5.10], [344.16,5.18], [347.55,5.14], [354.17,5.22] ]) # 5 March\n",
    "\n",
    "data_IT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 2.62976154e+02, -1.02172672e+05, -1.58510009e+02, -1.25290121e+08,\n",
       "       -1.58831083e+02])"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "popt2, pcov2 = scopt.curve_fit(f, t2, data_IT[:,0], p0=guess2)\n",
    "popt2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEGCAYAAACevtWaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXyU1dn/8c8FBBJkiQWMLEpQRKUgAaJSrMpmQVxwFxWrfVTqUlArPG617hXFasXHalFsXYsWISCi/GSJiAjIJqto2CQBlcVggknIcn5/3BMayDZJZjLb9/165eXMvc11nHBxOPe5r2POOUREJLo0CHUAIiISeEruIiJRSMldRCQKKbmLiEQhJXcRkSjUKNQBALRu3dolJyfX6tz9+/dzxBFHBDagMKc2xwa1OTbUpc3Lly/f7ZxrU9G+sEjuycnJLFu2rFbnpqen069fv8AGFObU5tigNseGurTZzLZVtk/DMiIiUUjJXUQkCim5i4hEobAYc69IYWEhmZmZ5OfnV3lcy5Yt2bBhQz1FFR4qanN8fDwdOnQgLi4uRFGJSDgJ2+SemZlJ8+bNSU5OxswqPS4nJ4fmzZvXY2Shd3ibnXPs2bOHzMxMOnXqFMLIRCRchO2wTH5+Pq1ataoysYvHzGjVqlW1/8oRkdgRtskdUGKvAf2/EpGywjq5i4hEreIiWPA0zXIygnJ5JXcRkfq2dzP8ayjMe5Q2uxYF5SOiJrmnrczijHHz6HTPB5wxbh5pK7MCct3vvvuO4cOHc/zxx9O1a1eGDh3K119/HZBrVyYtLY1HHnnk4Ov169cf3DdmzBg++eSToH6+iASJc7D8NXjx1/DDV3DpJLYc99ugfFS1yd3M4s1sqZl9aWbrzOxh33Yzs8fN7Gsz22Bmo8tsn2BmGWa22sx6BSXyMtJWZnHv1DVkZefhgKzsPO6duqbOCd45x8UXX0y/fv3YtGkT69ev5y9/+Qvff/99YAKvxFNPPcWtt94KlE/uo0aN4plnngnq54tIEOTugslXw/ujoUNvuHURdL8saB/nT8+9ABjgnOsBpABDzKwPcD1wDHCSc+5kYLLv+HOBE3w/I4EXAx304cbP3kheYfEh2/IKixk/e2Odrjt//nzi4uK4+eabD25LSUmhZ8+eDBw4kF69etG9e3emT58OeAWAzjvvPHr06EG3bt145513AFi+fDlnn302vXv3ZvDgwezcuROACRMm0LVrV0455RSGDx8OwNdff02TJk1o3bo1ixYtYsaMGYwdO5aUlBQ2bdpEx44d+fHHH/nuu+/q1DYRqUdfzYK/94GMuTD4Cbh2OrTsENSPrHaeu/MWWc31vY3z/TjgFuBq51yJ77gffMcMA173nbfYzBLNrK1zbmfAo/fZkZ1Xo+3+Wrt2Lb179y63PT4+nmnTptGiRQt2795Nnz59uPDCC/noo49o164dH3zwAQD79u2jsLCQUaNGMX36dNq0acM777zD/fffz6uvvsq4cePYsmULTZo0ITs7G4DPPvuMXr28f+z07duXCy+8kPPPP5/LLvvv3/A9evTgs88+49JLL61T+0QkyApyYfa9sOJ1OLo7XDITjjq5Xj7ar4eYzKwhsBzoDLzgnFtiZscDV5rZxcAuYLRz7hugPbC9zOmZvm07D7vmSLyePUlJSaSnpx/ymS1btiQnJ6fa2IqLizm6RRN2/lRQbt/RLZr4dY3K5Ofnc+DAgXLXKCws5J577mHRokU0aNCArKwsNm3aRKdOnfj444+58847GTJkCH379mXFihWsXbuWgQMHHow3KSmJnJwcunbtypVXXsl5553H+eefT8OGDdm6dSstWrQ4+JmFhYXk5eUdEkOrVq3YvHlzubjy8/PL/X+MFrm5uVHbtsqozZGtxb4NnLzhb8Tnf8+3x17K1uSrcOu/h/WHDusGq81+JXfnXDGQYmaJwDQz6wY0AfKdc6lmdgnwKnAmUNGEa1fBNScCEwFSU1Pd4SUvN2zY4NeTpzk5Odx97sncO3XNIUMzCXENufvck+v09Grv3r2ZOXNmuWv861//Yt++faxcuZK4uDiSk5Np1KgRvXr1YsWKFcyaNYtHH32U3/zmN1x88cX88pe/5PPPPy93/dmzZ7NgwQJmzJjB008/zbp160hMTGTfvn0HPzMuLo6EhIRDYjhw4ABHHnlkubji4+Pp2bNnrdsbzlQKNjZERZuLDsAnT8KqZ7yhl6s+pGPHX9GxksOD1eYazZZxzmUD6cAQvB75e75d04BTfK8z8cbiS3UAdtQpympc1LM9T1zSnfaJCRjQPjGBJy7pzkU929fpugMGDKCgoICXX3754LYvvviCbdu2cdRRRxEXF8f8+fPZts0rqbxjxw6aNm3KiBEjGDNmDCtWrODEE09k165dB5N7YWEh69ato6SkhO3bt9O/f3+eeuopsrOzyc3N5eSTTyYj47/zXps3b16uh56RkUG3bt3q1DYRCYJdG2HSIPj0aUi5Gm7+DDr+KiShVNtzN7M2QKFzLtvMEoBBwJNAGjAAr8d+NlA6P3AG8AczmwycDuwL5nh7qYt6tq9zMj+cmTFt2jTuuOMOxo0bR3x8PMnJyTz00EOMHj2a1NRUUlJSOOmkkwBYs2YNY8eOpUGDBsTFxfHiiy/SuHFjpkyZwujRo9m3bx9FRUXccccddOnShREjRrBv3z6cc9x5550kJiZy1llncdddd+Gcw8wYPnw4N910ExMmTGDKlCkce+yxbN68mdTU1IC2VUTqoKQEvngZPv4zND4CrnwTTr4gpCH5MyzTFnjNN+7eAHjXOTfTzBYCb5nZnXg3XG/0HT8LGApkAD8Dvwt82PWnXbt2vPvuu+W2VzTMkpyczODBg8ttT0lJYcGCBeW2L1y4sNy2pk2bMmjQIObOncugQYM444wzDpkKOW3aNIYNG0ajRmFb800ktvy0A9Juhc3z4YTBcOHz0Dwp1FH5NVtmNVBuINc3RHNeBdsdcFtAootR9913H0uWLKlwX1FREaNGjarniESkQmvfg5l/hOIDcP7foPf1ECZ1ntT9C0NJSUlceOGFFe67/PLL6zQDSEQCIC8bZo2FNe9C+1S4ZCK0Oj7UUR1CyV1EpCY2fwJpt0DOd9D/fvj1H6Fh+KXS8ItIRCQcFebD3Edg8QvQ6gS48WNoX/4hR3+lrcxi/OyNZGXn0X7xPMYOPjGgk0KU3EVEqrNzNUwdCbs2wKk3wTmPQOOmtb5caT2s0mdzSuthAQFL8FFTFTIYGjZsSEpKysGfrVu3smzZMkaPHg14Dx8sWhSccp0iEgZKimHhs/DyAMj7Ea55D857uk6JHYJXD6ss9dyrkJCQwKpVqw7ZlpycfHCOeXp6Os2aNaNv376hCE9EgmnPJm+K4/bF0HWYNxum6S8Cculg1cMqSz33GkpPT+f8889n69atvPTSSzz77LOkpKTw6aefhjo0EQmEkhJY+jK89GtvGObiiXD5awFL7ADtEhNqtL02IqPn/uE98N2aCnclFBfV7k710d3h3HFVHpKXl0dKSgoAnTp1Ytq0aQf3JScnc/PNN9OsWTPGjBlT888XkfCTvR2m3wZbPoHOg7wHklq0C/jHjB18YoX1sMYOPjFgnxEZyT1EKhqWEZEo5Bysegs+uhdcCVzwHPS6LmgPJJXeND04WyYxIUZny1TRw87LyalT5UcRiXE538H7t8PXH0HymTDsBTiyshqOgVNaDytYVSEjI7mHqebNm/PTTz+FOgwRqa2178EHd0FhHgx5Ek4bCQ2i41ZkdLQiRC644AKmTZumG6oikWb/Hnj3OpjyP9CqM9y8EPrcHDWJHdRzr1Jubm65bf369Tv4T6guXbqwevXqeo5KROrkq1neItV52TDwQeg7OizLB9RV9LVIRKQiedneDdMv3/Zmy12bBkdH76I3Su4iEv0y5sKMUd7N07P+F84aC40aV3lKae2XHdl5tAvCbJZgC+vkXroakVTPK6MvIocoyIWPH4Blr0LrE/0u9lUftV+CLWzvHsTHx7Nnzx4lLT8459izZw/x8fGhDkUkfGz9DF46A5b9E/qOgt8v8LuKY33Ufgm2sO25d+jQgczMTHbt2lXlcfn5+TGX1Cpqc3x8PB06dAhRRCJhpDAP5j0Gn/vmq//uwxovUl0ftV+CLWyTe1xcHJ06dar2uPT0dHr2LLcKYFSLxTaL+CVzOaTdDLu/hlNvhEEPQ5NmNb5Mu8QEsipI5IGs/RJsYTssIyLit6IDMPdRmHQOHNgP106D8/5aq8QOXu2XhLiGh2wLdO2XYAvbnruISHXSVmYx7cOPuDv/Obo22Ma2Yy+m49XPQXzLOl23bO0XzZYREalH05dvYfv0x3jFppFtzbjxwF18tuU0ntiQy0U965bc4b+1XyKVhmVEJPLsXE3XmRcxqsEUPig5nd8UPMmckt4RN6MlmNRzF5HIUXQAFoyHhc/QsqQZNxX+kY9LUg85JJJmtASTkruIRIYdKyHtNvhhHfS4it9+dS5fFZRPYZE0oyWYNCwjIuGtqADmPgIvD4S8vXDVO3DxS9w8JDXiZ7QEk3ruIhK+spZ7vfVdGyDlGhj8OCQcCUTHjJZgUnIXkfBTmA/pT8CiCdDsaLhmCpxwTrnDIn1GSzBVOyxjZvFmttTMvjSzdWb28GH7nzez3DLvm5jZO2aWYWZLzCw58GGLSNTa/gX840z47G/QcwTctrjCxC5V86fnXgAMcM7lmlkcsNDMPnTOLTazVCDxsONvAH50znU2s+HAk8CVgQ1bRKJOYR7HbfonfDIDmreDEe9B50GhjipiVdtzd57Snnmc78eZWUNgPPC/h50yDHjN93oKMNBUt1dEqvLtEnjp1xy7PQ16XQe3fq7EXkfmT0ldXyJfDnQGXnDO3W1mtwMNnHPPmlmuc66Z79i1wBDnXKbv/SbgdOfc7sOuORIYCZCUlNR78uTJtWpAbm4uzZrVrn5EpFKbY0MstLlBcQGdtrxJh8z3KWjShlUdbyC/XZ9Qh1Wv6vI99+/ff7lzLrWifX7dUHXOFQMpZpYITDOzs4DLgX4VHF5RL73c3yDOuYnARIDU1FRXui5pTaWnp1PbcyOV2hwboqXNla5otG0RTL8T9m6GU28kftBD5H++PCraXBPB+p5rNFvGOZdtZulAf7xefIZvxKWpmWU45zoDmcAxQKaZNQJaAnsDGrWIRISKVjR6ZOoXdF/zF47f/BYkHgvXvQ+dzgpxpNGn2uRuZm2AQl9iTwAGAU86544uc0yuL7EDzACuAz4HLgPmOS2nJBKTDl/R6HTbwFP2Dzpu/gFO+z0M/HOty/JK1fzpubcFXvONuzcA3nXOzazi+EnAG2aWgddjH173MEUkWIK5EHRpnZem5HN3o39zXaOP2VqSxJUFD/DO0DEB+QypWLXJ3Tm3Gqhy2Z/Sm6m+1/l44/EiEuaCvRB0u8QEkn9ayrhGr9DedjOp6FzGF11Bq8TDZ1BLoOkJVZEYVtVC0HVO7nk/8vZRr9MxfxqbStpyxYEHWOZOUv2XeqLkLhLDgrYQ9PrpMGssHffvZuMJN3Hzt4PYeqCY9qr/Um+U3EViWMAXgs75DmaNgQ3vw9GnwDX/4cS2PZhfxzil5lTyVySGBWwhaOdgxRvwwmnwzccw6CG4aT607RGwWKVm1HMXiWEBKZu7dwu8fzts+QQ6ngEXTIDWnas/T4JKyV0kxtW6bG5JMSx+EeY9Bg0awfnPQq/roYEGBMKBkruI1Nz362DGKG8xjS5D4LxnoKVukoYTJXcR8V9RASx4GhY+A/GJcOkk6HYpqPBr2FFyFxH/fLvE663v3ginXAmDn4AjWoU6KqmEkruIVK0g11ugeulEaNG+0iXvJLwouYtI5TLmwPt3wr7tcNpNvkJfzUMdlfhByV0kAgSzuFeFft4Ls++DL/8NrbvA/3wEx8bWIhqRTsldJMwFu7jXIZyDddPgw/+FvB/hzDFw1liIiw/s50jQaUKqSJirqrhXQP20AyZfDVN+542tj0yHgQ8osUco9dxFwlzQinuVKimBFa/Bx3+G4gNwzqPQ51ZoqPQQyfTtiYS5gBf3KuuHr7zSAdsXQ/KZcMFz0Or4ul9XQk7DMiJhLmDFvcoqzId5j8NLv/bmrQ97wVvLVIk9aqjnLhLmAlLcq6wtn8LMO2BPBnS/Agb/BZq1CWDEEg6U3EUiQK2Le5X18174fw/AqjchsSOMmAqdBwYmQAk7Su4i0c45WPMf+Oheb3rjGXfA2XdD46ahjkyCSMldJJrt3QIf/BE2zYP2veG3aXB091BHJfVAyV0kGhUXwucvQPo4aNAQzh0Pp97gvZaYoOQuEm0yl8P7o+H7tXDS+XDuU6q1HoOU3EWiRUEOzH3Uq97Y/Gi48k04+YJQRyUhouQuEg2++gBmjfVKCJx2Ewx4AOJbhDoqCSEld5FI9tMOL6l/NROO+iVc/hocc2qoo5IwoOQuEolKSmDZJJjzMJQUwsAHoe8oaBgX6sgkTFSb3M0sHlgANPEdP8U596CZvQWkAoXAUuD3zrlCMzPgOWAo8DNwvXNuRbAaIBJzvl/n1YPJ/AKO6wfnPwu/OC7UUUmY8ae2TAEwwDnXA0gBhphZH+At4CSgO5AA3Og7/lzgBN/PSODFQActEpMK87ye+j/Ogr2b4eKJcG2aErtUqNqeu3POAbm+t3G+H+ecm1V6jJktBTr43g4DXvedt9jMEs2srXNuZ2BDF4khGXPggzHw4xZIucYry6vFqaUKflWFNLOGZrYK+AH42Dm3pMy+OOBa4CPfpvbA9jKnZ/q2iUhN/bQT/nM9vHmp9wDSb2fARX9XYpdq+XVD1TlXDKSYWSIwzcy6OefW+nb/HVjgnPvU994qusThG8xsJN6wDUlJSaSnp9c0dgByc3NrfW6kUptjgCum9eZpFH06hQYlRWxLvppvj70E962Db9NDHV3QxNz3TPDaXKPZMs65bDNLB4YAa83sQaAN8Psyh2UCx5R53wHYUcG1JgITAVJTU12/fv1qFHip9PR0antupFKbo1zWcph5J+z8Eo4fAEOfplOr4+kU6rjqQUx9zz7BanO1wzJm1sbXY8fMEoBBwFdmdiMwGLjKOVdS5pQZwG/N0wfYp/F2iXZpK7M4Y9w8Ot3zAWeMm0fayqyaXyQvGz64C14eCDnfs67rWK8srxbQkFrwp+feFnjNzBri/WXwrnNuppkVAduAz73Zj0x1zj0CzMKbBpmBNxXyd0GJXCRMpK3M4t6paw4uYp2Vnce9U9cA+FeD3TlYMwVm3wc/74bTfw/972fX4hVgFY1yilTPn9kyq4GeFWyv8FzfLJnb6h6aSGQYP3vjwcReKq+wmPGzN1af3HdneCV5t3wC7XrBNf+BdilBjFZihZ5QFamjHRUsXl3VdsBbw3ThM7DwWWiUAOf9FXr/TiV5JWCU3EXqqF1iAlkVJPJ2iQkVn1B2znr3K+A3j0HzpCBHKbHGr3nuIlK5sYNPJCHu0B53QlxDxg4+8dADy81Znw6XvqzELkGhnrtIHZWOq4+fvZEd2Xm0S0xg7OAT/zveXlIMS1+GeY9B8QHofz+ccTs0ahLCqCXaKbmLBMBFPdtXfPP0kDnrA2HoeE1tlHqh5C4SDHnZMO9R+GISNEuCy/8FXS/S1EapN0ruIoFUbs76zdD/Pq2KJPVOyV0kUHZ/4z1hqjnrEgaU3EXq6sB+WDAeFv0fxDXVnHUJC0ruIrXlHKyfDrPvh58yvTrrgx6CZkeFOjIRJXeRWtn9jbcw9eb5kNQdLpsEx/YJdVQiBym5i9TE4UMw546H1P+BhvqjJOFFv5ESM9JWZlX+oFF1nIMNM+Cj+zQEIxFByV1iQp3K8moIRiKQkrvEhFqV5T2wHxY8DYue1xCMRBz9lkpMqFFZ3sOHYHpcDec8rCEYiShK7hIT/C7LuzsDPhwLm+ZpCEYimkr+Skyotizvgf0w52H4ex/IXAbnPgUj05XYJWKp5y4xodKyvCntYP0M+OheDcFIVFFyl5hRrizv7gx48xLfEEw3uPQV6Pir0AUoEkBK7hJ7DuyHT/8Kn02AuARvCCb1Bs2Ckaii32aJHeVmwVwF5zyiIRiJSkruEht+2AAf3u2V49UQjMQAJXeJbnnZkD4Olk6EJs1h6NNeOV4NwUiU02+4RKeSYlj5Bsx9BH7eC6m/g/5/giNahToykXqh5C5ho06Fvcr6don3INLOL+HYX8G5T0LbHoEPWCSMKblLWKhTYa9SP+2EOQ/C6negeTu4dBJ0u1SLUktMUnKXsFCrwl6ligpg8d/hk/FQUghnjoEz/wiNjwhixCLhrdrkbmbxwAKgie/4Kc65B82sEzAZ+AWwArjWOXfAzJoArwO9gT3Alc65rUGKX6JEjQp7lfX1bPjoHti7GU48DwY/Br84LggRikQWf2rLFAADnHM9gBRgiJn1AZ4EnnXOnQD8CNzgO/4G4EfnXGfgWd9xIlUqV8Crmu3szoC3Loe3rwBrCCPeg6veVmIX8ak2uTtPru9tnO/HAQOAKb7trwEX+V4P873Ht3+gmQY9pWrVFvYqVZADH//ZK/C17XP4zeNwyyLoPKgeoxUJf36NuZtZQ2A50Bl4AdgEZDvninyHZAKlA6Ptge0AzrkiM9sHtAJ2BzBuiTKVFvYqHW93JbDq394N09zvoecIGPigni4VqYQ55/w/2CwRmAb8Gfinb+gFMzsGmOWc625m64DBzrlM375NwGnOuT2HXWskMBIgKSmp9+TJk2vVgNzcXJo1a1arcyNVrLW5+U/fcNzGlzhyfwY/Ne/CNyfcRE6LLqEOK+hi7XsGtbmm+vfvv9w5l1rRvhrNlnHOZZtZOtAHSDSzRr7eewdgh++wTOAYINPMGgEtgb0VXGsiMBEgNTXV9evXryahHJSenk5tz41UMdPm3F0w92FY+SYH4lrCRS/S4pTh9G4QG8sQxMz3XIbaHDjV/ikxsza+HjtmlgAMAjYA84HLfIddB0z3vZ7he49v/zxXk38eiBQXwud/h+d7w5f/hr5/YMnpL0LK1RAjiV2krvzpubcFXvONuzcA3nXOzTSz9cBkM3sMWAlM8h0/CXjDzDLweuzDgxC3RKuMuTD7Ptj1FRw/EIaMgzZdKE5PD3VkIhGl2uTunFsN9Kxg+2bgtAq25wOXByQ6iR27v4HZ98M3s+HITnDVZOgyRE+XitSSnlCV0Mr7ET55yqvaGNcUznkUTv89NGoS6shEIpqSu4RGcREs/yfMfxzy90Gv66D//dCsTagjE4kKSu5S/8qOqyef6Y2rH90t1FGJRBUld6k/h4+rD38bThyqcXWRIFByl+DTuLpIvVNyl+DRuLpIyCi5S3BoXF0kpJTcJbA0ri4SFpTcpUYqXedU4+oiYUXJXfxW0Tqnf5q6ik6b36ZHxgsaVxcJI0ru4rfD1zk9s8FqHrA36LImS+PqImFGJfbEb6XrmR5nO5gUN543Go+jMUWMPPBHuO59JXaRMKKeu/ita8sDXLH/La5pOJefacLjhVfzWvFg2iS20A1TkTCj5C7VK8yHJS+SVvw01vBn3i4eyHNFl7CHlhWvcyoiIafkLpUrKYF1U2HOw7DvW+K6DGFOh9v4x6Ii9mbn0f7wdU5FJGwouUvFtn3uPYS0YwUc3R2GzYDjzmYQMOisUAcnItVRcpdD7dkEcx6EDe9D83Zw0YtwynAtbycSYZTcxfPzXu8hpC9egYaNof+f4Fe3QeOmoY5MRGpByT3WFRXA0pdhwVNQkAM9r/UeQmqeFOrIRKQOlNxjlXOwPg3mPAQ/boXOg7ySAUldQx2ZiASAknss2r7UK+6VuRSO+iWMmAqdB4Y6KhEJICX3WPLjVq+nvm4aNEuCC5+HlGugQcNQRyYiAabkHgvyfoRP/wpL/gHWEM6+G/qOhibNQh2ZiASJkns0KzoAy16FT8ZBXrbXSx9wP7RoF+rIRCTIlNyjUNqKTJZ8+DojC16jU4Pv+KF1H4767Xhoe0qoQxOReqInU6LMgrkzOWb6JTxR+CRFNOT6A2M5+/s7SPuuVahDE5F6pJ57tNj9Dcx5iLO+mskPJHJv4Q28W9yPYhpCSQnjZ29UDRiRGKLkHulyvoP0cbDidYhL4K+Fl/NK8bnkEX/IYaW12EUkNlQ7LGNmx5jZfDPbYGbrzOx23/YUM1tsZqvMbJmZnebbbmY2wcwyzGy1mfUKdiNiUkEOzHscJvSElW/AqTfA6FVMbX51ucQO0C4xIQRBikio+DPmXgTc5Zw7GegD3GZmXYGngIedcynAn33vAc4FTvD9jAReDHjUsay4kHZZH8BzKV7JgC6D4balMHQ8NGvD2MEnkhB36Lx11VwXiT3VDss453YCO32vc8xsA9AecEAL32EtgR2+18OA151zDlhsZolm1tZ3Hamt0nIBcx+hy97N3pql5zwM7XsfcljpuPr42RvZkZ1HO9VcF4lJ5uVgPw82SwYWAN3wEvxswPD+BdDXObfNzGYC45xzC33nzAXuds4tO+xaI/F69iQlJfWePHlyrRqQm5tLs2bR/TBOy+x1HL/pX7TI+Zr9TY9lXbsr+bn9GTG1tF0sfM+HU5tjQ13a3L9//+XOudQKdzrn/PoBmgHLgUt87ycAl/peXwHM8b3+APh1mfPmAr2runbv3r1dbc2fP7/W54a979c799YVzj3YwrmnT3JuxRvOFRdFd5sroTbHBrW5ZoBlrpK86tdsGTOLA94D3nLOTfVtvg643ff6P8ArvteZwDFlTu/Af4dsxB/7siD9L7DqbWjcHAY+CH1ugTjdFBUR/1Sb3M3MgEnABufcM2V27QDOBtKBAcA3vu0zgD+Y2WTgdGCf03i7f/L3wcJnYfGL4Erg9FvgrDHQ9BehjkxEIow/PfczgGuBNWa2yrftPuAm4DkzawTk4xs/B2YBQ4EM4GfgdwGNOBoVFcAXk2DBeMjbC92vgAF/giM7hjoyEYlQ/syWWYh307QivQ/f4BsHuq2OccWGkhJY+x7MewSyv4Xj+nszYNr2CHVkIhLh9IRqqGya7y1EvfNLOLq7FswQkYBScq9vmcth7kOwZQG0PBYungjdLx9A870AAAnnSURBVIcGquEmIoGjjFJfdm2Ed0bAKwMoyFrLs41u4MTvH+WMD9uQ9qXuN4tIYKnnHmzZ273FMla9DXFHsOGkPzBiXSp7ChsDkJWdx71T1wDoKVIRCRj13INl/x5vEerne8Pqd71pjbev4satAw8m9lJ5hcWMn70xRIGKSDRSzz3QCnLg87/DouehcD+kXA1n3wOJ3nNdlZXeVUleEQkkJfdAKSqAZf/05qr/vBtOvgAGPABtDq3G2C4xgawKErlK8opIIGlYpq5KimHVv+H5VPjobjjqZLhxLlz5ZrnEDqgkr4jUC/Xca8s52DgL5j4KuzZA2xS48DnvQaQqqjWqJK+I1Acl99rYuhDmPASZX0CrznD5a9B1mN8leC/q2V7JXESCSsm9JnZ+CXMfgYw50LwdXDABUq6BhvrfKCLhRVnJH3s2wbzHYN1USDgSznkUTrtJJXhFJGwpuVflp53wyZOw4nVo1ATOGgt9R0F8y1BHJiJSJSX3iuzfAwufgS9e8WbDnHoDnDkGmieFOjIREb8ouZeVlw2f/5+3WEbhz3DKlXD23fCLTqGOTESkRpTcAQpyYclLsGiCtxrSLy+GfvdWOE9dRCQSxHZyL8yHZZPg02e8p0q7DIH+90PbU0IdmYhIncRmci86ACvfgAVPQ84OOK4f9P8THHNqqCMTEQmI2EruJcWw+h1IHwfZ2+CY0+GSidDpzFBHJiISULGR3EtKYH0apD8Bu7/21ig976/QeZDfT5WKiESS6E7uzsHXH8G8x+H7NdDmJLjiDa9iYzVJPW1lluq/iEjEis7k7hxsTveeKs1aBkd2gktehm6XQoOG1Z6etjKLe6euIa+wGNBqSSISeaIvuX+72KvUuG0htOjgq/9yNTSM8/sS42dvPJjYS5WulqTkLiKRIHqS+46V3vBLxsdwxFFw7lPQ+3qvbEBNL6XVkkQkwkV8cj8idxu8MwI2vO8V9Rr0sFfUq/ERtb6mVksSkUgX2SsxrX6X1GW3w6Z074nS27+EX99Rp8QOWi1JRCJfZPfcjx/At8deSsfh46HpLwJ2Wa2WJCKRLrKT+xGt2XLctXQMYGIvpdWSRCSSVTssY2bHmNl8M9tgZuvM7PYy+0aZ2Ubf9qfKbL/XzDJ8+wYHK3gREamYPz33IuAu59wKM2sOLDezj4EkYBhwinOuwMyOAjCzrsBw4JdAO2COmXVxzhVXcn0REQmwanvuzrmdzrkVvtc5wAagPXALMM45V+Db94PvlGHAZOdcgXNuC5ABnBaM4EVEpGLmnPP/YLNkYAHQzfff6cAQIB8Y45z7wsz+D1jsnHvTd84k4EPn3JTDrjUSGAmQlJTUe/LkybVqQG5uLs2aNavVuZFKbY4NanNsqEub+/fvv9w5l1rRPr9vqJpZM+A94A7n3E9m1gg4EugDnAq8a2bHARUVbSn3N4hzbiIwESA1NdX169fP31AOkZ6eTm3PjVRqc2xQm2NDsNrs1zx3M4vDS+xvOeem+jZnAlOdZylQArT2bT+mzOkdgB2BC1lERKrjz2wZAyYBG5xzz5TZlQYM8B3TBWgM7AZmAMPNrImZdQJOAJYGOnAREamcP8MyZwDXAmvMbJVv233Aq8CrZrYWOABc57wB/HVm9i6wHm+mzW2aKSMiUr+qTe7OuYVUPI4OMKKScx4HHq9DXCIiUgcR+4Rq6WIaWdl5tF88T+UBRETKiMjkrsU0RESqFpFVIataTENERCI0uWsxDRGRqkVkcq9s0QwtpiEi4onI5K7FNEREqhaRN1TLLqaRlZ1Hey2mISJyiIhM7vDfxTRisRaFiEh1InJYRkREqqbkLiIShZTcRUSikJK7iEgUUnIXEYlCNVpmL2hBmO0CttXy9NZ4deRjidocG9Tm2FCXNnd0zrWpaEdYJPe6MLNlla0hGK3U5tigNseGYLVZwzIiIlFIyV1EJApFQ3KfGOoAQkBtjg1qc2wISpsjfsxdRETKi4aeu4iIHEbJXUQkCkVMcjezIWa20cwyzOyeCvY3MbN3fPuXmFly/UcZWH60+Y9mtt7MVpvZXDPrGIo4A6m6Npc57jIzc2YW8dPm/GmzmV3h+67Xmdnb9R1joPnxu32smc03s5W+3++hoYgzUMzsVTP7wczWVrLfzGyC7//HajPrVecPdc6F/Q/QENgEHAc0Br4Euh52zK3AS77Xw4F3Qh13PbS5P9DU9/qWWGiz77jmwAJgMZAa6rjr4Xs+AVgJHOl7f1So466HNk8EbvG97gpsDXXcdWzzWUAvYG0l+4cCHwIG9AGW1PUzI6XnfhqQ4Zzb7Jw7AEwGhh12zDDgNd/rKcBAM7N6jDHQqm2zc26+c+5n39vFQId6jjHQ/PmeAR4FngLy6zO4IPGnzTcBLzjnfgRwzv1QzzEGmj9tdkAL3+uWwI56jC/gnHMLgL1VHDIMeN15FgOJZta2Lp8ZKcm9PbC9zPtM37YKj3HOFQH7gFb1El1w+NPmsm7A+5s/klXbZjPrCRzjnJtZn4EFkT/fcxegi5l9ZmaLzWxIvUUXHP60+SFghJllArOAUfUTWsjU9M97tSJlJaaKeuCHz+H055hI4nd7zGwEkAqcHdSIgq/KNptZA+BZ4Pr6Cqge+PM9N8IbmumH96+zT82sm3MuO8ixBYs/bb4K+Jdz7q9m9ivgDV+bS4IfXkgEPH9FSs89EzimzPsOlP9n2sFjzKwR3j/lqvpnULjzp82Y2SDgfuBC51xBPcUWLNW1uTnQDUg3s614Y5MzIvymqr+/29Odc4XOuS3ARrxkH6n8afMNwLsAzrnPgXi8AlvRyq8/7zURKcn9C+AEM+tkZo3xbpjOOOyYGcB1vteXAfOc705FhKq2zb4hin/gJfZIH4eFatrsnNvnnGvtnEt2ziXj3We40Dm3LDThBoQ/v9tpeDfPMbPWeMM0m+s1ysDyp83fAgMBzOxkvOS+q16jrF8zgN/6Zs30AfY553bW6Yqhvotcg7vNQ4Gv8e6y3+/b9gjeH27wvvz/ABnAUuC4UMdcD22eA3wPrPL9zAh1zMFu82HHphPhs2X8/J4NeAZYD6wBhoc65npoc1fgM7yZNKuA34Q65jq299/ATqAQr5d+A3AzcHOZ7/gF3/+PNYH4vVb5ARGRKBQpwzIiIlIDSu4iIlFIyV1EJAopuYuIRCEldxGRKKTkLlIJM0s0s1tDHYdIbSi5i1QuEa/aqEjEUXIXqdw44HgzW2Vm40MdjEhN6CEmkUr4FnyZ6ZzrFuJQRGpMPXcRkSik5C4iEoWU3EUql4NXZlgk4ii5i1TCObcH+MzM1uqGqkQa3VAVEYlC6rmLiEQhJXcRkSik5C4iEoWU3EVEopCSu4hIFFJyFxGJQkruIiJR6P8DKMBJM7T3yw8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.plot(t2, data_IT[:, 0], 'o',  label='Cases(t)')\n",
    "# plt.plot(t2, data_IT[:, 1], 'o', label='Deaths(t)')\n",
    "plt.plot(t2, f(t2, *popt2), '-',  label='Fit')\n",
    "\n",
    "plt.legend(loc='best')\n",
    "plt.xlabel('t')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}