Compartmental_model.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Compartmental_model.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AJX1jgLnX70s"
      },
      "source": [
        "# Epidemiology model\n",
        "\n",
        "https://nbviewer.jupyter.org/github/pyro-ppl/pyro/blob/sir-tutorial-ii/tutorial/source/epi_regional.ipynb?fbclid=IwAR3Gv8tLuiEjOmZh7-NQUa_ggm_QUqtSc5TxRZ0_pSxVA7Y3lWWzSFGKjrA \n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "z6UoAzRe1pMh"
      },
      "source": [
        "!git clone https://github.com/pyro-ppl/pyro.git"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ryMxWMvbD8Nc"
      },
      "source": [
        "%cd /content/pyro\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "w8MT-jR48mLX"
      },
      "source": [
        "!pip install .[extras]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7AUI1jmXcX4u"
      },
      "source": [
        "import os\n",
        "import logging\n",
        "import urllib.request\n",
        "from collections import OrderedDict\n",
        "\n",
        "import pandas as pd\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "import torch\n",
        "import pyro\n",
        "import pyro.distributions as dist\n",
        "from pyro.contrib.epidemiology import CompartmentalModel, binomial_dist, infection_dist\n",
        "from pyro.ops.tensor_utils import convolve\n",
        "\n",
        "%matplotlib inline\n",
        "pyro.enable_validation(True)           \n",
        "torch.set_default_dtype(torch.double) \n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hzpnb36feNgS"
      },
      "source": [
        "  ## Model without Policies\n",
        "  "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bsaJcjx6xVLo"
      },
      "source": [
        "class CovidModel(CompartmentalModel):\n",
        "    def __init__(self, population, new_cases, new_recovered, new_deaths):\n",
        "        '''\n",
        "        population (int) – Total population = S + E + I + R.\n",
        "        '''\n",
        "        assert len(new_cases) == len(new_recovered) == len(new_deaths)\n",
        "\n",
        "        compartments = (\"S\", \"E\", \"I\", \"D\")  # R is implicit.\n",
        "        duration = len(new_cases)\n",
        "        super().__init__(compartments, duration, population)\n",
        "\n",
        "        self.new_cases = new_cases\n",
        "        self.new_deaths = new_deaths\n",
        "        self.new_recovered = new_recovered\n",
        "        \n",
        "\n",
        "    def global_model(self):\n",
        "        tau_i = pyro.sample(\"rec_time\", dist.Normal(15.0, 3.0))\n",
        "        tau_e = pyro.sample(\"incub_time\", dist.Normal(5.0, 1.0))\n",
        "        # R0 = pyro.sample(\"R0\", dist.LogNormal(0., 1.))\n",
        "        R0 = pyro.sample(\"R0\", dist.Normal(2.5, 0.5))\n",
        "        rho = pyro.sample(\"rho\", dist.Beta(10, 10))  # About 50% response rate.\n",
        "        mort_rate = pyro.sample(\"mort_rate\", dist.Beta(2, 50))  # About 2% mortality rate.\n",
        "        rec_rate = pyro.sample(\"rec_rate\",dist.Beta(10, 10)) # About 50% recovery rate.\n",
        "        return R0, tau_e, tau_i, rho, mort_rate, rec_rate\n",
        "\n",
        "    def initialize(self, params):\n",
        "        # Start with a single infection.\n",
        "        return {\"S\": self.population - 1, \"E\": 0, \"I\": 1, \"D\": 0}\n",
        "\n",
        "    def transition(self, params, state, t):\n",
        "        R0, tau_e, tau_i, rho, mort_rate, rec_rate = params\n",
        "\n",
        "        # Sample flows between compartments.\n",
        "        S2E = pyro.sample(\"S2E_{}\".format(t),\n",
        "                            infection_dist(individual_rate=R0 / tau_i,\n",
        "                                          num_susceptible=state[\"S\"],\n",
        "                                          num_infectious=state[\"I\"],\n",
        "                                          population=self.population))\n",
        "        E2I = pyro.sample(\"E2I_{}\".format(t),\n",
        "                            binomial_dist(state[\"E\"], 1 / tau_e )) \n",
        "        I2R = pyro.sample(\"I2R_{}\".format(t),\n",
        "                            binomial_dist(state[\"I\"], 1 / tau_i))\n",
        "        I2D = pyro.sample(\"I2D_{}\".format(t),\n",
        "                          binomial_dist(state[\"I\"], mort_rate / tau_i))\n",
        "\n",
        "        # Update compartments with flows.\n",
        "        state[\"S\"] = state[\"S\"] - S2E \n",
        "        state[\"E\"] = state[\"E\"] + S2E - E2I\n",
        "        state[\"I\"] = state[\"I\"] + E2I - I2R - I2D\n",
        "        state[\"D\"] = state[\"D\"] + I2D\n",
        "\n",
        "        # Condition on observations.\n",
        "        t_is_observed = isinstance(t, slice) or t < self.duration\n",
        "        pyro.sample(\"new_cases_{}\".format(t),\n",
        "                    binomial_dist(S2E, rho),\n",
        "                    obs=self.new_cases[t] if t_is_observed else None)\n",
        "        pyro.sample(\"new_deaths_{}\".format(t),\n",
        "                      binomial_dist(I2D, 1),\n",
        "                      obs=self.new_deaths[t] if t_is_observed else None)\n",
        "        pyro.sample(\"new_recovered_{}\".format(t),\n",
        "                    binomial_dist(I2R, rho),\n",
        "                    obs=self.new_recovered[t] if t_is_observed else None)\n",
        "        \n",
        "    def compute_flows(self, prev, curr, t):\n",
        "        S2E = prev[\"S\"] - curr[\"S\"]  # S can only go to E.\n",
        "        I2D = curr[\"D\"] - prev[\"D\"]  # D can only have come from I.\n",
        "        # We deduce the remaining flows by conservation of mass:\n",
        "        #   curr - prev = inflows - outflows\n",
        "        E2I = prev[\"E\"] - curr[\"E\"] + S2E\n",
        "        I2R = prev[\"I\"] - curr[\"I\"] + E2I - I2D\n",
        "        return {\n",
        "            \"S2E_{}\".format(t): S2E,\n",
        "            \"E2I_{}\".format(t): E2I,\n",
        "            \"I2D_{}\".format(t): I2D,\n",
        "            \"I2R_{}\".format(t): I2R,\n",
        "        }"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kgKCZvRfMi_3"
      },
      "source": [
        "## Create Country"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "koX5yGHrsuib"
      },
      "source": [
        "# function to make the time series of confirmed and daily confirmed cases for a specific country\n",
        "def create_country (country, start_date, end_date, state = False) : \n",
        "\n",
        "    url = 'https://raw.githubusercontent.com/assemzh/ProbProg-COVID-19/master/full_grouped.csv'\n",
        "    data = pd.read_csv(url)\n",
        "\n",
        "    data.Date = pd.to_datetime(data.Date)\n",
        "\n",
        "    if state :\n",
        "        df = data.loc[data[\"Province/State\"] == country, [\"Province/State\", \"Date\", \"Confirmed\", \"Deaths\", \"Recovered\", \"Active\", \"New cases\", \"New deaths\", \"New recovered\"]]\n",
        "    else : \n",
        "        df = data.loc[data[\"Country/Region\"] == country, [\"Country/Region\", \"Date\", \"Confirmed\", \"Deaths\", \"Recovered\", \"Active\", \"New cases\", \"New deaths\", \"New recovered\"]]\n",
        "    df.columns = [\"country\", \"date\", \"confirmed\", \"deaths\", \"recovered\", \"active\", \"new_cases\", \"new_deaths\", \"new_recovered\"]\n",
        "\n",
        "    # group by country and date\n",
        "    df = df.groupby(['country','date'])['confirmed', 'deaths', 'recovered',\"active\", \"new_cases\", \"new_deaths\", \"new_recovered\"].sum().reset_index()\n",
        "\n",
        "    # convert date string to datetime\n",
        "    df.date = pd.to_datetime(df.date)\n",
        "    df = df.sort_values(by = \"date\")\n",
        "    df = df[df.date >= start_date]\n",
        "    df = df[df.date <= end_date]\n",
        "\n",
        "    active = df['active'].tolist()\n",
        "    recovered = df['recovered'].tolist()\n",
        "    deaths = df['deaths'].tolist()\n",
        "    new_cases = df['new_cases'].tolist()\n",
        "    new_recovered = df['new_recovered'].tolist()\n",
        "    new_deaths = df['new_deaths'].tolist()\n",
        "    \n",
        "    active = torch.tensor(list(map(float, active))).view(len(active),1) \n",
        "    recovered = torch.tensor(list(map(float, recovered))).view(len(recovered),1) \n",
        "    deaths = torch.tensor(list(map(float, deaths))).view(len(deaths),1) \n",
        "    new_cases = torch.tensor(list(map(float, new_cases))).view(len(new_cases),1)  \n",
        "    new_recovered = torch.tensor(list(map(float, new_recovered))).view(len(new_recovered),1) \n",
        "    new_deaths = torch.tensor(list(map(float, new_deaths))).view(len(new_deaths),1) \n",
        "\n",
        "\n",
        "    return_data = {\n",
        "        'active':active,\n",
        "        'recovered':recovered,\n",
        "        'deaths':deaths,\n",
        "        'new_cases':new_cases,\n",
        "        'new_recovered': new_recovered,\n",
        "        'new_deaths':new_deaths }\n",
        "        \n",
        "    return return_data\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nxRxbNr8zt3O"
      },
      "source": [
        "## Get data for countries\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "MzNNysenCiWz"
      },
      "source": [
        "# Parameters\n",
        "country = \"Japan\"\n",
        "start_date = \"2020-02-01\" \n",
        "end_date = \"2020-04-01\"\n",
        "population = 126500000\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "V4O6XUbM9Ff3"
      },
      "source": [
        "data = create_country(country, start_date, end_date)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XTUPDWD9e_9o"
      },
      "source": [
        "##Train the model using MCMC.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "bX4Jq5qmC6Ke",
        "outputId": "dfce3586-740f-44e2-9324-25406d50cf8f"
      },
      "source": [
        "%%time\n",
        "model = CovidModel(population, data[\"new_cases\"], data[\"new_recovered\"], data[\"new_deaths\"] )\n",
        "mcmc = model.fit_mcmc(num_samples=500, warmup_steps = 200)\n",
        "mcmc.summary()"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "INFO \t Running inference...\n",
            "Warmup:   0%|          | 0/700 [00:00, ?it/s]INFO \t Heuristic init: R0=1.96, incub_time=4.5, mort_rate=0.00865, rec_rate=0.55, rec_time=16.7, rho=0.176\n",
            "Sample: 100%|██████████| 700/700 [03:19,  3.51it/s, step size=1.19e-03, acc. prob=0.919]\n"
          ],
          "name": "stderr"
        },
        {
          "output_type": "stream",
          "text": [
            "\n",
            "                     mean       std    median      5.0%     95.0%     n_eff     r_hat\n",
            "             R0      1.92      0.01      1.92      1.91      1.93      2.54      2.55\n",
            " auxiliary[0,0] 126499843.01      0.34 126499842.95 126499842.51 126499843.61      4.48      1.14\n",
            " auxiliary[0,1] 126499707.58      0.59 126499707.44 126499706.71 126499708.60      3.22      1.83\n",
            " auxiliary[0,2] 126499567.86      3.73 126499567.41 126499562.87 126499573.83      2.58      2.41\n",
            " auxiliary[0,3] 126499427.98      6.65 126499426.72 126499417.67 126499438.37      2.63      2.40\n",
            " auxiliary[0,4] 126499285.36     10.77 126499282.20 126499270.45 126499303.10      2.78      2.18\n",
            " auxiliary[0,5] 126499140.60     16.45 126499134.80 126499119.30 126499168.53      2.75      2.21\n",
            " auxiliary[0,6] 126498994.65     21.93 126498988.06 126498965.60 126499030.59      2.76      2.22\n",
            " auxiliary[0,7] 126498849.08     27.55 126498841.25 126498812.46 126498893.46      2.75      2.24\n",
            " auxiliary[0,8] 126498703.46     32.95 126498694.98 126498657.16 126498755.75      2.77      2.21\n",
            " auxiliary[0,9] 126498555.53     38.07 126498547.62 126498503.32 126498616.82      2.76      2.24\n",
            "auxiliary[0,10] 126498406.90     42.98 126498398.01 126498345.57 126498477.19      2.78      2.21\n",
            "auxiliary[0,11] 126498258.83     47.41 126498250.89 126498191.94 126498336.96      2.77      2.23\n",
            "auxiliary[0,12] 126498109.82     52.12 126498101.67 126498037.05 126498193.50      2.76      2.24\n",
            "auxiliary[0,13] 126497961.07     56.21 126497951.79 126497881.99 126498050.66      2.77      2.22\n",
            "auxiliary[0,14] 126497811.85     60.11 126497801.63 126497726.97 126497905.09      2.76      2.25\n",
            "auxiliary[0,15] 126497662.11     63.79 126497653.57 126497572.73 126497762.75      2.74      2.27\n",
            "auxiliary[0,16] 126497511.32     67.68 126497503.48 126497415.50 126497617.71      2.72      2.30\n",
            "auxiliary[0,17] 126497361.15     71.21 126497351.45 126497261.15 126497471.38      2.70      2.32\n",
            "auxiliary[0,18] 126497209.38     74.65 126497198.74 126497103.91 126497325.11      2.69      2.33\n",
            "auxiliary[0,19] 126497056.62     78.13 126497046.60 126496945.46 126497176.56      2.66      2.38\n",
            "auxiliary[0,20] 126496904.74     81.08 126496894.31 126496786.90 126497026.69      2.65      2.40\n",
            "auxiliary[0,21] 126496752.75     83.58 126496742.24 126496629.41 126496878.76      2.65      2.38\n",
            "auxiliary[0,22] 126496599.67     86.43 126496590.91 126496471.54 126496729.43      2.64      2.40\n",
            "auxiliary[0,23] 126496444.75     88.94 126496437.84 126496311.34 126496578.31      2.63      2.41\n",
            "auxiliary[0,24] 126496287.59     91.25 126496278.95 126496151.63 126496425.63      2.62      2.42\n",
            "auxiliary[0,25] 126496129.19     93.76 126496122.69 126495986.95 126496270.25      2.60      2.45\n",
            "auxiliary[0,26] 126495969.68     95.78 126495963.18 126495824.39 126496114.69      2.60      2.45\n",
            "auxiliary[0,27] 126495810.39     97.00 126495804.13 126495659.87 126495956.57      2.60      2.45\n",
            "auxiliary[0,28] 126495648.05     99.42 126495642.88 126495500.61 126495803.57      2.59      2.47\n",
            "auxiliary[0,29] 126495482.17    102.52 126495476.77 126495327.64 126495640.66      2.57      2.50\n",
            "auxiliary[0,30] 126495314.55    105.78 126495311.68 126495155.44 126495475.63      2.56      2.52\n",
            "auxiliary[0,31] 126495145.72    107.74 126495141.51 126494978.34 126495308.38      2.55      2.53\n",
            "auxiliary[0,32] 126494973.03    108.77 126494968.30 126494806.64 126495139.86      2.55      2.55\n",
            "auxiliary[0,33] 126494796.17    109.95 126494793.20 126494625.08 126494964.19      2.56      2.53\n",
            "auxiliary[0,34] 126494613.90    111.12 126494612.29 126494441.68 126494783.64      2.55      2.55\n",
            "auxiliary[0,35] 126494428.69    111.60 126494428.62 126494254.55 126494599.00      2.54      2.56\n",
            "auxiliary[0,36] 126494241.16    111.63 126494239.10 126494064.60 126494410.23      2.55      2.54\n",
            "auxiliary[0,37] 126494049.79    111.77 126494048.04 126493875.01 126494217.59      2.55      2.54\n",
            "auxiliary[0,38] 126493853.19    113.47 126493850.78 126493678.49 126494023.04      2.53      2.58\n",
            "auxiliary[0,39] 126493653.51    115.06 126493654.87 126493476.20 126493825.01      2.52      2.61\n",
            "auxiliary[0,40] 126493452.24    117.10 126493451.79 126493269.80 126493626.60      2.52      2.61\n",
            "auxiliary[0,41] 126493247.33    118.28 126493244.29 126493063.30 126493422.59      2.51      2.63\n",
            "auxiliary[0,42] 126493036.42    119.32 126493034.72 126492849.65 126493213.87      2.50      2.66\n",
            "auxiliary[0,43] 126492823.12    119.55 126492820.03 126492638.50 126493000.20      2.50      2.67\n",
            "auxiliary[0,44] 126492603.06    121.96 126492604.03 126492408.43 126492782.45      2.51      2.67\n",
            "auxiliary[0,45] 126492377.22    123.96 126492377.46 126492203.34 126492578.74      2.51      2.63\n",
            "auxiliary[0,46] 126492145.11    124.15 126492143.26 126491967.79 126492346.29      2.52      2.61\n",
            "auxiliary[0,47] 126491906.38    124.38 126491903.36 126491731.95 126492104.03      2.52      2.63\n",
            "auxiliary[0,48] 126491660.08    122.90 126491656.55 126491486.79 126491857.68      2.53      2.60\n",
            "auxiliary[0,49] 126491408.14    121.43 126491405.26 126491236.00 126491604.28      2.53      2.59\n",
            "auxiliary[0,50] 126491154.32    119.57 126491158.17 126490988.48 126491349.37      2.51      2.64\n",
            "auxiliary[0,51] 126490902.57    117.02 126490905.74 126490738.36 126491089.58      2.51      2.63\n",
            "auxiliary[0,52] 126490647.59    114.14 126490650.82 126490482.12 126490826.25      2.50      2.64\n",
            "auxiliary[0,53] 126490389.80    109.80 126490396.98 126490235.21 126490560.56      2.50      2.65\n",
            "auxiliary[0,54] 126490131.51    105.23 126490139.63 126489981.08 126490292.88      2.50      2.63\n",
            "auxiliary[0,55] 126489871.99    100.48 126489882.12 126489730.22 126490024.63      2.50      2.62\n",
            "auxiliary[0,56] 126489613.40     97.03 126489621.68 126489471.72 126489756.17      2.50      2.62\n",
            "auxiliary[0,57] 126489356.12     95.69 126489361.52 126489214.53 126489498.29      2.51      2.60\n",
            "auxiliary[0,58] 126489106.11     96.65 126489110.87 126488965.04 126489247.90      2.52      2.58\n",
            "auxiliary[0,59] 126488859.47     98.07 126488866.70 126488714.73 126489003.86      2.53      2.56\n",
            "auxiliary[0,60] 126488616.68    100.16 126488626.56 126488470.18 126488765.94      2.56      2.51\n",
            " auxiliary[1,0]    143.44      0.45    143.65    142.59    143.93      3.41      1.61\n",
            " auxiliary[1,1]    230.55      0.22    230.56    230.15    230.85     13.33      1.01\n",
            " auxiliary[1,2]    280.31      0.54    280.34    279.43    281.13      8.41      1.07\n",
            " auxiliary[1,3]    313.28      2.05    313.09    309.73    316.28      3.64      1.69\n",
            " auxiliary[1,4]    336.19      3.01    335.42    332.08    341.36      2.81      2.11\n",
            " auxiliary[1,5]    358.54      0.91    358.49    357.15    360.06      5.32      1.45\n",
            " auxiliary[1,6]    359.90      1.35    359.71    358.19    362.37      4.03      1.55\n",
            " auxiliary[1,7]    377.75      0.56    377.70    377.05    378.83      6.56      1.00\n",
            " auxiliary[1,8]    384.13      1.13    384.27    382.54    385.94      3.88      1.03\n",
            " auxiliary[1,9]    410.96      2.34    410.33    407.67    414.33      2.66      3.25\n",
            "auxiliary[1,10]    425.11      1.04    424.95    423.57    426.53      2.88      2.29\n",
            "auxiliary[1,11]    431.11      0.56    431.29    430.09    431.76      4.06      1.84\n",
            "auxiliary[1,12]    442.09      0.76    442.07    441.06    443.23      3.95      1.54\n",
            "auxiliary[1,13]    455.47      0.69    455.53    454.71    456.50      5.86      1.09\n",
            "auxiliary[1,14]    456.20      0.49    456.19    455.49    457.14     10.48      1.04\n",
            "auxiliary[1,15]    457.87      0.96    458.16    456.41    459.19      3.58      1.75\n",
            "auxiliary[1,16]    466.77      1.20    466.50    464.98    468.76      5.26      1.20\n",
            "auxiliary[1,17]    450.11      0.49    450.02    449.45    451.07      7.35      1.04\n",
            "auxiliary[1,18]    460.99      0.52    461.01    460.22    461.90      9.96      1.00\n",
            "auxiliary[1,19]    490.70      0.55    490.53    489.95    491.60      7.39      1.40\n",
            "auxiliary[1,20]    499.72      0.80    499.86    498.23    500.83      4.68      1.00\n",
            "auxiliary[1,21]    491.06      0.71    490.85    490.07    492.09      6.75      1.09\n",
            "auxiliary[1,22]    496.38      1.30    496.89    493.72    497.81      3.16      1.80\n",
            "auxiliary[1,23]    528.14      0.87    528.39    526.65    529.43      4.80      1.33\n",
            "auxiliary[1,24]    541.36      0.92    541.62    539.77    542.55      3.65      1.33\n",
            "auxiliary[1,25]    564.22      0.84    564.15    562.93    565.67      4.51      1.76\n",
            "auxiliary[1,26]    603.03      1.09    602.84    601.36    605.01      7.25      1.28\n",
            "auxiliary[1,27]    614.56      0.43    614.59    613.95    615.21      9.68      1.29\n",
            "auxiliary[1,28]    630.75      1.66    630.28    628.69    633.64      2.82      2.10\n",
            "auxiliary[1,29]    669.88      0.68    669.74    668.85    671.08      6.52      1.36\n",
            "auxiliary[1,30]    685.86      1.81    686.20    683.18    688.22      2.60      3.05\n",
            "auxiliary[1,31]    694.05      1.03    694.01    692.67    695.87      4.42      1.65\n",
            "auxiliary[1,32]    728.61      3.84    730.30    722.44    733.03      2.50      2.71\n",
            "auxiliary[1,33]    767.86      3.50    769.56    763.08    771.91      2.50      3.10\n",
            "auxiliary[1,34]    823.21      1.17    823.65    821.04    824.55      3.76      1.67\n",
            "auxiliary[1,35]    861.48      1.81    861.76    858.76    864.45      3.44      1.99\n",
            "auxiliary[1,36]    883.27      1.10    883.52    881.18    884.84      6.00      1.02\n",
            "auxiliary[1,37]    926.68      1.31    926.66    924.64    928.68      3.40      2.34\n",
            "auxiliary[1,38]   1015.41      2.92   1015.12   1011.11   1019.96      3.02      1.95\n",
            "auxiliary[1,39]   1085.89      3.69   1084.85   1080.83   1091.17      2.83      2.22\n",
            "auxiliary[1,40]   1101.16      1.62   1101.06   1098.93   1103.96      3.76      1.77\n",
            "auxiliary[1,41]   1105.15      2.80   1104.07   1102.41   1111.20      5.62      1.16\n",
            "auxiliary[1,42]   1162.25      1.92   1161.40   1160.28   1165.75      3.29      1.82\n",
            "auxiliary[1,43]   1197.62      1.98   1197.44   1194.33   1200.68      7.88      1.03\n",
            "auxiliary[1,44]   1246.34      1.80   1246.11   1243.90   1249.35      2.94      2.22\n",
            "auxiliary[1,45]   1288.54      1.63   1288.56   1286.00   1290.88      6.32      1.31\n",
            "auxiliary[1,46]   1313.38      2.00   1313.00   1310.56   1316.83      5.31      1.25\n",
            "auxiliary[1,47]   1339.84      2.37   1339.08   1337.66   1345.38      4.81      1.30\n",
            "auxiliary[1,48]   1412.79      1.33   1412.82   1410.61   1415.09     19.06      1.01\n",
            "auxiliary[1,49]   1506.06      3.27   1506.49   1501.62   1510.87      4.50      1.01\n",
            "auxiliary[1,50]   1601.76      1.75   1601.72   1599.21   1604.94      5.41      1.42\n",
            "auxiliary[1,51]   1603.07      3.43   1602.00   1599.41   1609.33      4.08      1.24\n",
            "auxiliary[1,52]   1665.80      1.74   1665.73   1663.84   1669.07      8.94      1.26\n",
            "auxiliary[1,53]   1752.78      5.80   1749.69   1746.55   1761.69      3.01      1.84\n",
            "auxiliary[1,54]   1820.87      4.41   1820.60   1814.80   1828.09      3.56      1.32\n",
            "auxiliary[1,55]   1901.04      6.92   1903.07   1890.60   1910.54      2.94      2.12\n",
            "auxiliary[1,56]   1995.35      7.24   1994.23   1986.53   2010.15      3.01      1.98\n",
            "auxiliary[1,57]   2089.19      9.92   2088.38   2072.20   2105.10      3.05      1.97\n",
            "auxiliary[1,58]   2112.25      8.82   2112.89   2100.37   2125.74      2.51      2.86\n",
            "auxiliary[1,59]   2139.01      8.64   2135.22   2128.22   2150.59      2.38      3.60\n",
            "auxiliary[1,60]   2097.98     32.56   2100.35   2044.57   2145.76      2.58      2.51\n",
            " auxiliary[2,0]      0.68      0.00      0.68      0.67      0.68      5.68      1.04\n",
            " auxiliary[2,1]     33.38      0.03     33.38     33.32     33.42      9.43      1.16\n",
            " auxiliary[2,2]     87.18      0.25     87.10     86.86     87.55      2.64      2.66\n",
            " auxiliary[2,3]    147.93      0.33    147.93    147.37    148.39      3.63      1.68\n",
            " auxiliary[2,4]    215.20      0.44    215.13    214.55    215.88      6.72      1.00\n",
            " auxiliary[2,5]    280.64      1.58    280.83    277.92    283.03      2.81      2.23\n",
            " auxiliary[2,6]    364.38      1.08    364.45    362.26    365.81      4.56      1.45\n",
            " auxiliary[2,7]    424.96      1.52    424.86    422.23    427.03      5.95      1.34\n",
            " auxiliary[2,8]    493.39      0.63    493.41    492.43    494.34     13.54      1.33\n",
            " auxiliary[2,9]    541.27      2.30    541.72    537.86    544.49      2.56      2.46\n",
            "auxiliary[2,10]    600.58      1.81    601.22    597.18    602.61      3.63      1.56\n",
            "auxiliary[2,11]    663.72      1.62    663.31    661.90    666.94      3.65      1.65\n",
            "auxiliary[2,12]    720.89      1.47    720.69    718.70    723.70      8.07      1.22\n",
            "auxiliary[2,13]    771.45      2.21    772.10    767.94    774.50      3.04      2.25\n",
            "auxiliary[2,14]    834.64      2.52    834.92    830.92    838.49      3.46      1.95\n",
            "auxiliary[2,15]    895.67      1.76    895.67    893.21    898.80      3.56      2.08\n",
            "auxiliary[2,16]    948.67      2.28    947.84    945.80    952.21      3.51      1.66\n",
            "auxiliary[2,17]   1024.21      1.27   1024.14   1022.30   1026.12      6.68      1.13\n",
            "auxiliary[2,18]   1073.34      1.46   1073.29   1070.97   1075.97      8.67      1.11\n",
            "auxiliary[2,19]   1103.16      2.31   1103.17   1099.57   1106.84      5.67      2.04\n",
            "auxiliary[2,20]   1148.69      2.49   1149.06   1144.95   1152.19      3.12      2.89\n",
            "auxiliary[2,21]   1211.69      1.47   1211.90   1208.97   1213.84      8.05      1.02\n",
            "auxiliary[2,22]   1261.75      2.86   1261.30   1257.68   1266.19      3.65      1.84\n",
            "auxiliary[2,23]   1286.69      0.97   1286.75   1284.92   1288.02      9.99      1.00\n",
            "auxiliary[2,24]   1338.50      1.98   1338.75   1335.06   1341.60      5.18      1.06\n",
            "auxiliary[2,25]   1381.41      1.96   1381.41   1377.94   1384.39      4.65      1.40\n",
            "auxiliary[2,26]   1408.21      2.92   1409.05   1403.70   1412.09      4.18      1.40\n",
            "auxiliary[2,27]   1462.08      2.61   1462.50   1457.22   1465.29      3.70      1.68\n",
            "auxiliary[2,28]   1519.71      2.60   1519.77   1515.24   1524.13      5.51      1.36\n",
            "auxiliary[2,29]   1562.06      4.33   1562.56   1556.03   1568.38      2.63      2.57\n",
            "auxiliary[2,30]   1632.14      6.88   1630.89   1622.25   1641.72      2.48      2.49\n",
            "auxiliary[2,31]   1712.80      8.52   1713.07   1701.59   1726.20      2.58      2.42\n",
            "auxiliary[2,32]   1773.66      3.99   1773.77   1767.46   1779.67      2.80      2.61\n",
            "auxiliary[2,33]   1837.47      4.29   1836.77   1832.48   1845.51      4.81      1.36\n",
            "auxiliary[2,34]   1894.21      5.61   1893.00   1886.13   1903.44      2.92      2.21\n",
            "auxiliary[2,35]   1972.06      6.25   1971.16   1964.24   1982.13      2.65      2.82\n",
            "auxiliary[2,36]   2069.76      3.30   2070.15   2063.21   2073.62      4.67      1.47\n",
            "auxiliary[2,37]   2150.89      1.89   2150.60   2147.74   2153.84      7.01      1.21\n",
            "auxiliary[2,38]   2192.80      1.91   2192.86   2189.38   2195.47     10.17      1.18\n",
            "auxiliary[2,39]   2256.97      6.01   2256.35   2248.50   2266.01      2.50      2.75\n",
            "auxiliary[2,40]   2378.46      5.18   2378.31   2371.63   2386.39      2.58      2.57\n",
            "auxiliary[2,41]   2515.58      3.60   2514.54   2510.88   2521.72      4.37      1.58\n",
            "auxiliary[2,42]   2606.25      2.98   2606.86   2601.06   2610.32      4.02      2.08\n",
            "auxiliary[2,43]   2720.84      3.86   2721.09   2715.01   2727.21      2.80      2.34\n",
            "auxiliary[2,44]   2829.46      3.57   2829.68   2823.58   2834.85      5.16      1.80\n",
            "auxiliary[2,45]   2950.61      4.88   2950.34   2942.43   2958.40      3.41      1.81\n",
            "auxiliary[2,46]   3095.99      4.00   3095.22   3090.43   3103.41      6.02      1.30\n",
            "auxiliary[2,47]   3247.07      4.91   3247.15   3239.49   3255.24      4.77      1.76\n",
            "auxiliary[2,48]   3359.73      3.04   3359.49   3354.69   3364.76     10.41      1.13\n",
            "auxiliary[2,49]   3458.17      5.28   3456.92   3450.48   3466.76      3.83      1.48\n",
            "auxiliary[2,50]   3556.18      6.97   3556.26   3545.94   3567.85      3.43      1.74\n",
            "auxiliary[2,51]   3746.70     10.47   3748.13   3729.73   3761.57      3.10      2.06\n",
            "auxiliary[2,52]   3878.85     18.39   3878.78   3852.05   3911.78      2.88      2.16\n",
            "auxiliary[2,53]   3989.51     18.64   3992.31   3961.23   4018.97      2.81      2.22\n",
            "auxiliary[2,54]   4119.69     32.45   4120.08   4072.76   4168.54      2.58      2.54\n",
            "auxiliary[2,55]   4239.10     31.01   4229.12   4192.72   4289.14      2.67      2.41\n",
            "auxiliary[2,56]   4343.09     35.74   4336.56   4287.34   4393.87      2.58      2.56\n",
            "auxiliary[2,57]   4445.77     36.14   4444.89   4390.51   4500.71      2.53      2.67\n",
            "auxiliary[2,58]   4611.05     37.54   4609.48   4557.25   4671.46      2.62      2.44\n",
            "auxiliary[2,59]   4768.43     38.13   4767.37   4709.81   4830.26      2.67      2.33\n",
            "auxiliary[2,60]   4989.34     14.37   4983.25   4970.67   5012.81      2.74      2.20\n",
            " auxiliary[3,0]      0.48      0.00      0.48      0.48      0.49      5.81      1.35\n",
            " auxiliary[3,1]      0.50      0.00      0.50      0.50      0.51      3.99      1.70\n",
            " auxiliary[3,2]      0.50      0.00      0.50      0.50      0.51      3.36      2.09\n",
            " auxiliary[3,3]      0.50      0.00      0.50      0.50      0.50      5.87      1.04\n",
            " auxiliary[3,4]      0.50      0.00      0.50      0.50      0.50      3.46      1.65\n",
            " auxiliary[3,5]      0.49      0.00      0.49      0.49      0.50      6.11      1.51\n",
            " auxiliary[3,6]      0.50      0.00      0.49      0.49      0.50      2.91      2.19\n",
            " auxiliary[3,7]      0.50      0.00      0.50      0.50      0.51      2.77      2.41\n",
            " auxiliary[3,8]      0.51      0.00      0.51      0.51      0.52      4.51      1.38\n",
            " auxiliary[3,9]      0.48      0.00      0.48      0.48      0.49      3.41      2.32\n",
            "auxiliary[3,10]      0.50      0.00      0.50      0.49      0.50      2.82      2.03\n",
            "auxiliary[3,11]      0.50      0.00      0.50      0.50      0.51      2.71      2.38\n",
            "auxiliary[3,12]      0.79      0.00      0.79      0.78      0.79      3.00      2.18\n",
            "auxiliary[3,13]      1.15      0.00      1.15      1.15      1.16      3.14      2.16\n",
            "auxiliary[3,14]      1.27      0.00      1.26      1.26      1.27      6.91      1.15\n",
            "auxiliary[3,15]      1.16      0.01      1.16      1.15      1.17      2.62      2.74\n",
            "auxiliary[3,16]      1.78      0.00      1.78      1.77      1.79      3.25      2.03\n",
            "auxiliary[3,17]      2.41      0.01      2.41      2.39      2.42      2.77      2.34\n",
            "auxiliary[3,18]      2.66      0.01      2.66      2.65      2.67      5.72      1.34\n",
            "auxiliary[3,19]      2.63      0.01      2.63      2.62      2.64      3.93      1.40\n",
            "auxiliary[3,20]      2.88      0.01      2.88      2.87      2.89      5.67      1.26\n",
            "auxiliary[3,21]      3.00      0.02      3.00      2.98      3.03      2.92      2.20\n",
            "auxiliary[3,22]      2.90      0.00      2.90      2.89      2.91      9.60      1.05\n",
            "auxiliary[3,23]      3.14      0.00      3.14      3.14      3.15     12.92      1.13\n",
            "auxiliary[3,24]      3.52      0.03      3.52      3.47      3.57      3.00      2.08\n",
            "auxiliary[3,25]      4.50      0.03      4.50      4.45      4.55      2.93      2.17\n",
            "auxiliary[3,26]      5.71      0.03      5.71      5.65      5.75      4.11      2.61\n",
            "auxiliary[3,27]      6.90      0.02      6.91      6.86      6.94      3.31      1.96\n",
            "auxiliary[3,28]      7.73      0.04      7.73      7.67      7.81      3.55      1.82\n",
            "auxiliary[3,29]      8.28      0.09      8.28      8.16      8.41      2.49      2.95\n",
            "auxiliary[3,30]      8.63      0.07      8.63      8.51      8.73      2.66      2.40\n",
            "auxiliary[3,31]      8.91      0.02      8.91      8.89      8.93      8.12      1.07\n",
            "auxiliary[3,32]      9.13      0.03      9.14      9.08      9.17      3.89      1.60\n",
            "auxiliary[3,33]      9.41      0.08      9.42      9.29      9.53      2.53      2.65\n",
            "auxiliary[3,34]      9.91      0.14      9.90      9.69     10.13      2.52      2.64\n",
            "auxiliary[3,35]     11.03      0.26     11.04     10.64     11.42      2.56      2.53\n",
            "auxiliary[3,36]     12.34      0.29     12.37     11.86     12.77      2.64      2.34\n",
            "auxiliary[3,37]     13.71      0.32     13.71     13.24     14.25      2.60      2.44\n",
            "auxiliary[3,38]     15.38      0.26     15.34     15.03     15.79      2.53      2.59\n",
            "auxiliary[3,39]     17.07      0.20     17.02     16.82     17.39      2.46      2.88\n",
            "auxiliary[3,40]     19.06      0.12     19.02     18.91     19.24      2.47      2.81\n",
            "auxiliary[3,41]     21.11      0.07     21.09     21.03     21.22      2.78      2.15\n",
            "auxiliary[3,42]     23.23      0.03     23.23     23.17     23.26      5.75      1.36\n",
            "auxiliary[3,43]     25.32      0.08     25.34     25.20     25.43      2.79      2.08\n",
            "auxiliary[3,44]     27.38      0.12     27.42     27.19     27.54      2.78      2.09\n",
            "auxiliary[3,45]     29.45      0.16     29.50     29.20     29.70      2.83      2.03\n",
            "auxiliary[3,46]     31.48      0.21     31.52     31.12     31.78      2.69      2.34\n",
            "auxiliary[3,47]     33.34      0.30     33.41     32.83     33.72      2.56      2.59\n",
            "auxiliary[3,48]     35.10      0.39     35.15     34.53     35.73      2.52      2.71\n",
            "auxiliary[3,49]     36.73      0.52     36.75     35.97     37.53      2.45      2.84\n",
            "auxiliary[3,50]     38.22      0.64     38.20     37.30     39.17      2.46      2.76\n",
            "auxiliary[3,51]     39.63      0.75     39.60     38.58     40.77      2.48      2.71\n",
            "auxiliary[3,52]     40.96      0.84     40.90     39.81     42.28      2.50      2.66\n",
            "auxiliary[3,53]     42.22      0.92     42.16     40.98     43.67      2.51      2.62\n",
            "auxiliary[3,54]     43.42      1.00     43.36     41.95     44.86      2.53      2.55\n",
            "auxiliary[3,55]     44.57      1.05     44.50     43.06     46.11      2.54      2.53\n",
            "auxiliary[3,56]     45.67      1.10     45.61     44.13     47.30      2.54      2.54\n",
            "auxiliary[3,57]     46.74      1.13     46.63     45.15     48.37      2.54      2.57\n",
            "auxiliary[3,58]     47.77      1.15     47.66     46.12     49.38      2.53      2.61\n",
            "auxiliary[3,59]     48.77      1.16     48.69     47.05     50.37      2.54      2.64\n",
            "auxiliary[3,60]     49.74      1.17     49.67     48.05     51.36      2.54      2.66\n",
            "     incub_time      4.52      0.00      4.52      4.52      4.53      3.31      1.72\n",
            "      mort_rate      0.01      0.00      0.01      0.01      0.01      2.71      2.24\n",
            "       rec_rate      0.55      0.00      0.55      0.55      0.55      9.73      1.00\n",
            "       rec_time     16.69      0.00     16.69     16.69     16.70      2.73      2.71\n",
            "            rho      0.26      0.00      0.26      0.26      0.27      2.58      2.48\n",
            "\n",
            "Number of divergences: 0\n",
            "CPU times: user 3min 18s, sys: 1.21 s, total: 3min 19s\n",
            "Wall time: 3min 19s\n"
          ],
          "name": "stdout"
        }
      ]
    }
  ]
}