.. _2020covidrst:

=======================
Algo - simulation COVID
=======================


.. only:: html

    **Links:** :download:`notebook <2020_covid.ipynb>`, :downloadlink:`html <2020_covid2html.html>`, :download:`python <2020_covid.py>`, :downloadlink:`slides <2020_covid.slides.html>`, :githublink:`GitHub|_doc/notebooks/td1a_home/2020_covid.ipynb|*`


Ou comment utiliser les mathématiques pour comprendre la propagation de
l’épidémie.

.. code:: ipython3

    from jyquickhelper import add_notebook_menu
    add_notebook_menu()






.. contents::
    :local:





.. code:: ipython3

    %matplotlib inline

Enoncé
------

On récupère les données du COVID par région et par âge et premier graphe
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

A cette adresse : `Données relatives à l’épidémie de COVID-19 en France
: vue
d’ensemble <https://www.data.gouv.fr/en/datasets/donnees-relatives-a-lepidemie-de-covid-19-en-france-vue-densemble/>`__

.. code:: ipython3

    from pandas import read_csv, to_datetime
    url = "https://www.data.gouv.fr/en/datasets/r/d3a98a30-893f-47f7-96c5-2f4bcaaa0d71"
    covid = read_csv(url, sep=",")
    covid['date'] = to_datetime(covid['date'])
    covid.tail()






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>date</th>
          <th>total_cas_confirmes</th>
          <th>total_deces_hopital</th>
          <th>total_deces_ehpad</th>
          <th>total_cas_confirmes_ehpad</th>
          <th>total_cas_possibles_ehpad</th>
          <th>patients_reanimation</th>
          <th>patients_hospitalises</th>
          <th>total_patients_gueris</th>
          <th>nouveaux_patients_hospitalises</th>
          <th>nouveaux_patients_reanimation</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>802</th>
          <td>2022-05-13</td>
          <td>NaN</td>
          <td>118431</td>
          <td>28826.0</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>1233.0</td>
          <td>19272.0</td>
          <td>653050.0</td>
          <td>756.0</td>
          <td>77.0</td>
        </tr>
        <tr>
          <th>803</th>
          <td>2022-05-14</td>
          <td>NaN</td>
          <td>118494</td>
          <td>28837.0</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>1213.0</td>
          <td>18901.0</td>
          <td>653734.0</td>
          <td>387.0</td>
          <td>42.0</td>
        </tr>
        <tr>
          <th>804</th>
          <td>2022-05-15</td>
          <td>NaN</td>
          <td>118508</td>
          <td>28844.0</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>1214.0</td>
          <td>18935.0</td>
          <td>653827.0</td>
          <td>145.0</td>
          <td>13.0</td>
        </tr>
        <tr>
          <th>805</th>
          <td>2022-05-16</td>
          <td>NaN</td>
          <td>118633</td>
          <td>28845.0</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>1199.0</td>
          <td>18742.0</td>
          <td>654775.0</td>
          <td>900.0</td>
          <td>94.0</td>
        </tr>
        <tr>
          <th>806</th>
          <td>2022-05-17</td>
          <td>NaN</td>
          <td>118723</td>
          <td>28845.0</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>1173.0</td>
          <td>18290.0</td>
          <td>655984.0</td>
          <td>883.0</td>
          <td>106.0</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    ax = covid.set_index("date").plot(
        title="Evolution des hospitalisations par jour", figsize=(14, 4))
    ax.set_yscale("log");



.. image:: 2020_covid_6_0.png


Il y a quelques valeurs manquantes même pour les séries aggrégées… Comme
je n’ai pas le courage de corriger les valeurs unes à unes, je prends un
autre fichier et quelques aberrations comme le nombre de décès qui
décroît, ce qui est rigoureusement impossible.

.. code:: ipython3

    from pandas import concat, to_datetime
    
    def extract_data(kind='deaths', country='France'):
        url = (
            "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/"
            "master/csse_covid_19_data/"
            "csse_covid_19_time_series/time_series_covid19_%s_global.csv" %
            kind)
        df = read_csv(url)
        eur = df[df['Country/Region'].isin([country])
                 & df['Province/State'].isna()]
        tf = eur.T.iloc[4:]
        tf.columns = [kind]
        return tf
    
    def extract_whole_data(kind=['deaths', 'confirmed', 'recovered'],
                           country='France'):
        population = {
            'France': 67e6,
        }
    
        total = population[country]
        dfs = []
        for k in kind:
            df = extract_data(k, country)
            dfs.append(df)
        conc = concat(dfs, axis=1)
        conc['infected'] = conc['confirmed'] - (conc['deaths'] + conc['recovered'])
        conc['safe'] = total - conc.drop('confirmed', axis=1).sum(axis=1)
        conc.index = to_datetime(conc.index)
        
        return conc
    
    covid = extract_whole_data()
    covid.tail()






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>deaths</th>
          <th>confirmed</th>
          <th>recovered</th>
          <th>infected</th>
          <th>safe</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>2022-11-29</th>
          <td>155309</td>
          <td>36762891</td>
          <td>0</td>
          <td>36607582</td>
          <td>30237109.0</td>
        </tr>
        <tr>
          <th>2022-11-30</th>
          <td>155384</td>
          <td>36829869</td>
          <td>0</td>
          <td>36674485</td>
          <td>30170131.0</td>
        </tr>
        <tr>
          <th>2022-12-01</th>
          <td>155459</td>
          <td>36897875</td>
          <td>0</td>
          <td>36742416</td>
          <td>30102125.0</td>
        </tr>
        <tr>
          <th>2022-12-02</th>
          <td>155526</td>
          <td>36955945</td>
          <td>0</td>
          <td>36800419</td>
          <td>30044055.0</td>
        </tr>
        <tr>
          <th>2022-12-03</th>
          <td>155526</td>
          <td>36955945</td>
          <td>0</td>
          <td>36800419</td>
          <td>30044055.0</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    import matplotlib.pyplot as plt
    fig, ax = plt.subplots(1, 2, figsize=(14, 4))
    covid[['confirmed', 'infected']].plot(title="Evolution de l'épidémie par jour", ax=ax[0])
    covid[['deaths', 'recovered']].plot(title="Evolution de l'épidémie par jour", ax=ax[1]);



.. image:: 2020_covid_9_0.png


Même aberration, un nombre décès qui décroît… Il faudrait comprendre
pourquoi pour savoir comment réparer les données. Ou on improvise. Pour
chaque observation :math:`X_t < X_{t-1}`, on calcule le ratio
:math:`\frac{X_{t}}{X_{t-1}}` et on multiplie toutes les observations
:math:`i<t` par ce ratio.

.. code:: ipython3

    from tqdm import tqdm  # pour avoir une barre de progression
    
    def correct_series(X):
        for t in range(1, X.shape[0]):
            if X[t-1] > 0 and X[t] == 0:
                X[t] = X[t-1]
                continue
            if X[t] >= X[t-1] and X[t] < X[t-1] + 200000:
                continue
            ratio = X[t] / X[t-1]
            for i in range(0, t):
                X[i] *= ratio
    
    covid_modified = covid.copy()
    for c in tqdm(covid.columns):
        values = covid_modified[c].values
        correct_series(values)
        covid_modified[c] = values
    
    covid_modified.tail()


.. parsed-literal::
    100%|██████████| 5/5 [00:00<00:00, 29.74it/s]






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>deaths</th>
          <th>confirmed</th>
          <th>recovered</th>
          <th>infected</th>
          <th>safe</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>2022-11-29</th>
          <td>155309</td>
          <td>36762891</td>
          <td>342253</td>
          <td>36607582</td>
          <td>30044055.0</td>
        </tr>
        <tr>
          <th>2022-11-30</th>
          <td>155384</td>
          <td>36829869</td>
          <td>342253</td>
          <td>36674485</td>
          <td>30044055.0</td>
        </tr>
        <tr>
          <th>2022-12-01</th>
          <td>155459</td>
          <td>36897875</td>
          <td>342253</td>
          <td>36742416</td>
          <td>30044055.0</td>
        </tr>
        <tr>
          <th>2022-12-02</th>
          <td>155526</td>
          <td>36955945</td>
          <td>342253</td>
          <td>36800419</td>
          <td>30044055.0</td>
        </tr>
        <tr>
          <th>2022-12-03</th>
          <td>155526</td>
          <td>36955945</td>
          <td>342253</td>
          <td>36800419</td>
          <td>30044055.0</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    fig, ax = plt.subplots(1, 2, figsize=(14, 4))
    covid_modified[['confirmed', 'infected']].plot(title="Evolution de l'épidémie par jour", ax=ax[0])
    covid_modified[['deaths', 'recovered']].plot(title="Evolution de l'épidémie par jour", ax=ax[1]);



.. image:: 2020_covid_12_0.png


C’est mieux.

.. code:: ipython3

    covid = covid_modified

On lisse.

.. code:: ipython3

    lisse = covid.rolling(7).mean()
    fig, ax = plt.subplots(1, 2, figsize=(14, 4))
    lisse[['confirmed', 'infected']].plot(title="Evolution de l'épidémie par jour", ax=ax[0])
    lisse[['deaths', 'recovered']].plot(title="Evolution de l'épidémie par jour", ax=ax[1]);



.. image:: 2020_covid_16_0.png


La série des *recovered* ne compte vraisemblablement que ceux qui sont
passés à l’hôpital. Il faudrait recouper avec d’autres données pour être
sûr. Ce sera pour un autre jour.

Modèle `SIRD <http://www.xavierdupre.fr/app/aftercovid/helpsphinx/model_covidsir.html>`__
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Pour en savoir plus `Modèles compartimentaux en
épidémiologie <https://fr.wikipedia.org/wiki/Mod%C3%A8les_compartimentaux_en_%C3%A9pid%C3%A9miologie>`__.
On classe la population en quatre catégories :


-  S : personnes non contaminées
-  I : nombre de personnes malades ou contaminantes
-  R : personnes guéries (recovered)
-  D : personnes décédées

Les gens changent de catégorie en fonction de l’évolution de l’épidémie
selon les équations qui suivent :


-  :math:`\frac{dS}{dt} = - \beta \frac{S I}{N}`
-  :math:`\frac{dI}{dt} = \frac{\beta S I}{N} - \mu I - \nu I`
-  :math:`\frac{dD}{dt} = \nu I`
-  :math:`\frac{dR}{dt} = \mu I`

:math:`\beta` est lié au taux de transmission, :math:`\frac{1}{\mu}` est
la durée moyenne jusqu’à guérison, :math:`\frac{1}{\nu}` est la durée
moyenne jusqu’au décès.

Q0 : une petite fonction pour dessiner
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Cette fonction servira à représenter graphiquement les résultats.

.. code:: ipython3

    from datetime import datetime, timedelta
    
    
    def plot_simulation(sim, day0=datetime(2020, 1, 1), safe=True,
                        ax=None, title=None, logy=False, two=False,
                        true_data=None):
        """
        On suppose que sim est une matrice (days, 4).
        :param sim: la simulation
        :param day0: le premier jour de la simulation (une observation par jour)
        :param safe: ajouter les personnes *safe* (non infectées), comme elles sont nombreuses,
            il vaut mieux aussi cocher *logy=True* pour que cela soit lisible
        :param ax: axes existant (utile pour superposer), None pour un créer un nouveau
        :param title: titre du graphe
        :param logy: échelle logarithmique sur l'axe des y
        :param two: faire deux graphes plutôt qu'un seul pour plus de visibilité
        :param true_data: vraies données à tracer également en plus de celle de la simulation
        :return: ax
        """
        df = DataFrame(sim, columns=['S', 'I', 'R', 'D'])
        # On ajoute des dates.
        df["date"] = [day0 + timedelta(d) for d in range(0, df.shape[0])]
        df = df.set_index("date")
        if true_data is None:
            tdf = None
        else:
            tdf = DataFrame(true_data, columns=['Sobs', 'Iobs', 'Robs', 'Dobs'])
            tdf["date"] = [day0 + timedelta(d) for d in range(0, tdf.shape[0])]
            tdf = tdf.set_index("date")
            
        if two:
            if ax is None:
                fig, ax = plt.subplots(1, 2, figsize=(14, 4))
            if safe:
                if tdf is not None:
                    tdf.drop(['Dobs'], axis=1).plot(ax=ax[0], logy=logy, linewidth=8)
                df.drop('D', axis=1).plot(ax=ax[0], title=title, logy=logy, linewidth=4)
            else:
                if tdf is not None:
                    tdf.drop(['Sobs', 'Dobs'], axis=1).plot(ax=ax[0], logy=logy, linewidth=8)
                df.drop(['S', 'D'], axis=1).plot(ax=ax[0], title=title, logy=logy, linewidth=4)
            if tdf is not None:
                tdf['Dobs'].plot(ax=ax[1], title=title, logy=logy, linewidth=8)
            df[['D']].plot(ax=ax[1], title='Décès', logy=logy, linewidth=4)
            ax[0].legend()
            ax[1].legend()
        else:
            if ax is None:
                fig, ax = plt.subplots(1, 1, figsize=(14, 4))
            if safe:
                if tdf is not None:
                    tdf.plot(ax=ax, title=title, logy=logy, linewidth=8)
                df.plot(ax=ax, title=title, logy=logy, linewidth=4)
            else:
                if tdf is not None:
                    tdf.drop(['Sobs'], axis=1).plot(ax=ax, title=title, logy=logy, linewidth=8)
                df.drop(['S'], axis=1).plot(ax=ax, title=title, logy=logy, linewidth=4)
            ax.legend()
        return ax

Q1 : écrire une fonction qui calcule la propagation
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

On suppose que :math:`\beta, \mu, \nu, S_0, I_0, R_0, D_0` sont connus.
On rappelle le modèle :


-  :math:`dS = - \beta \frac{S I}{N}`
-  :math:`dI = \frac{\beta S I}{N} - \mu I - \nu I`
-  :math:`dD = \nu I`
-  :math:`dR = \mu I`

.. code:: ipython3

    import numpy

.. code:: ipython3

    beta = 0.5
    mu = 1./14
    nu = 1./21
    S0 = 9990
    I0 = 10
    R0 = 0
    D0 = 0

Il faudra compléter le petit programme suivant :

.. code:: ipython3

    from pandas import DataFrame
    
    
    
    def simulation(beta, mu, nu, S0, I0, R0, D0, days=14):
        res = numpy.empty((days+1, 4), dtype=numpy.float64)
        res[0, :] = [S0, I0, R0, D0]
        N = sum(res[0, :])
        for t in range(1, res.shape[0]):
            dR = res[t-1, 1] * mu
            # ....
        return res
    
    
    sim = simulation(beta, mu, nu, S0, I0, R0, D0, 30)
    plot_simulation(sim);



.. image:: 2020_covid_25_0.png


Q2 : on veut estimer les paramètres du modèle, une fonction d’erreur ?
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

C’est compliqué parce que… les paramètres évoluent au cours du temps, en
fonction du comportement des gens, masque, pas masque, confinement,
reconfinement, température, manque de tests également… Tout d’abord les
vraies données.

.. code:: ipython3

    lisse_mars = lisse[30:]
    dates = lisse_mars.index
    france = numpy.zeros((lisse_mars.shape[0], 4), dtype=numpy.dtype)
    france[:, 3] = lisse_mars['deaths']
    france[:, 2] = lisse_mars['recovered']
    france[:, 0] = lisse_mars['safe']
    france[:, 1] = lisse_mars['infected']
    
    france[:5]




.. parsed-literal::
    array([[29963471.89949666, 1.1640875027069664, 3.626756310367156,
            0.9858981877075415],
           [29963471.89949666, 1.1640875027069664, 3.626756310367156,
            0.9858981877075415],
           [29963471.89949666, 1.1640875027069664, 3.626756310367156,
            0.9858981877075415],
           [29963471.89949666, 1.1640875027069664, 3.626756310367156,
            0.9858981877075415],
           [29963471.89949666, 1.1640875027069664, 4.533445387958945,
            0.9858981877075415]], dtype=object)



.. code:: ipython3

    plot_simulation(france, dates[0], safe=False, logy=True, title="Vraies données");



.. image:: 2020_covid_28_0.png


Et sur les derniers jours.

.. code:: ipython3

    plot_simulation(france[-60:], dates[-60], two=True, safe=False, title="Vraies données, derniers mois");



.. image:: 2020_covid_30_0.png


Bref, on part du principe que le modèle est plutôt fiable sur une courte
période de temps, on tire plein de paramètres aléatoires et on regarde
ce qui marche le mieux. Et pour comparer deux jeux de paramètres, il
faut donc une fonction d’erreur qu’on prendra comme égal à la somme des
erreurs de prédictions.

Maintenant il faut faire attention à ce qu’on compare. La simulation
calcule les catégories de population au temps *t*, mais pas toujours les
séries cumulées. La série des personnes contaminées est **transitoire
dans la simulation** et **cumulées dans les données récupérées**. La
première étape consiste à transformer les données simulées pour qu’elles
soient comparables aux données collectées.

.. code:: ipython3

    def simulation_cumulee(beta, mu, nu, S0, I0, R0, D0, days=14):
        # ...
        pass

Maintenant la fonction d’erreur :

.. code:: ipython3

    def error(data, simulation):
        # ... à compléter
        return 0

Q3 : optimisation
~~~~~~~~~~~~~~~~~

Pour optimiser, on tire des paramètres de façon aléatoire dans un
intervalle donné et on choisit ceux qui minimisent l’erreur.

.. code:: ipython3

    from tqdm import tqdm  # pour avoir une barre de progression
    
    
    def optimisation(true_data, i_range=(0, 0.2), beta_range=(0, 0.5),
                     mu_range=(0., 0.2), nu_range=(0., 0.2),
                     max_iter=1000, error_fct=error):
        N = sum(true_data[0, :])
        rnd = numpy.random.rand(max_iter, 4)
        for i, (a, b) in enumerate([i_range, beta_range, mu_range, nu_range]):
            rnd[:, i] = rnd[:, i] # à compléter ... 
            
        err_min = None
        for it in tqdm(range(max_iter)):
            i, beta, mu, nu = rnd[it, :]
            
            D0 = true_data[0, 3]
            # dI0 = 
            # S0 = 
            # I0 = 
            # R0 = 
            sim = simulation_cumulee(beta, mu, nu, S0, I0, R0, D0, days=true_data.shape[0] - 1)
            err = error_fct(true_data, sim)
            if err_min is None or err < err_min:
                # à compléter
                pass
        return best


Q4 : dessiner les résultats
~~~~~~~~~~~~~~~~~~~~~~~~~~~



Q5 : vérifier que cela marche sur des données synthétiques
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

On simule, on vérifie que l’optimisation retrouve les paramètres de la
simulation.


Q6 : sur des données réelles
~~~~~~~~~~~~~~~~~~~~~~~~~~~~


Réponses
--------

Q1 : propagation
~~~~~~~~~~~~~~~~

On l’applique aux données réelles.

.. code:: ipython3

    from datetime import datetime, timedelta
    from pandas import DataFrame
    
    
    def simulation(beta, mu, nu, S0, I0, R0, D0, days=14):
        res = numpy.empty((days+1, 4), dtype=numpy.float64)
        res[0, :] = [S0, I0, R0, D0]
        N = sum(res[0, :])
        for t in range(1, res.shape[0]):
            dR = res[t-1, 1] * mu
            dD = res[t-1, 1] * nu
            dI = res[t-1, 0] * res[t-1, 1] / N * beta
            res[t, 0] = res[t-1, 0] - dI
            res[t, 1] = res[t-1, 1] + dI - dR - dD
            res[t, 2] = res[t-1, 2] + dR
            res[t, 3] = res[t-1, 3] + dD
        return res    
    
    
    beta = 0.5
    mu = 1./14
    nu = 1./21
    S0 = 9990
    I0 = 10
    R0 = 0
    D0 = 0
    
    sim = simulation(beta, mu, nu, S0, I0, R0, D0, 60)
    plot_simulation(sim, dates[60], safe=False, two=True,
                    title="Simulation pour essayer");



.. image:: 2020_covid_47_0.png


Q2 : série cumulées et fonction erreur
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

On doit d’abord calculer la simulation qui modifie :math:`I_t` en
:math:`J_t` qui correspond à l’ensemble des personnes contaminées
jusqu’à présent.

.. code:: ipython3

    def simulation_cumulee(beta, mu, nu, S0, I0, R0, D0, days=14):
        res = numpy.empty((days+1, 4), dtype=numpy.float64)
        cum = numpy.empty((days+1, 1), dtype=numpy.float64)
        res[0, :] = [S0, I0, R0, D0]
        cum[0, 0] = I0
        N = sum(res[0, :])
        for t in range(1, res.shape[0]):
            dR = res[t-1, 1] * mu
            dD = res[t-1, 1] * nu
            dI = res[t-1, 0] * res[t-1, 1] / N * beta
            res[t, 0] = res[t-1, 0] - dI
            res[t, 1] = res[t-1, 1] + dI - dR - dD
            res[t, 2] = res[t-1, 2] + dR
            res[t, 3] = res[t-1, 3] + dD
            cum[t, 0] = cum[t-1, 0] + dI
        res[:, 1] = cum[:, 0]
        return res
    
    beta = 0.5
    mu = 1./14
    nu = 1./21
    S0 = 9990
    I0 = 10
    R0 = 0
    D0 = 0
    
    sim = simulation_cumulee(beta, mu, nu, S0, I0, R0, D0, 60)
    plot_simulation(sim, dates[60], safe=False, two=True,
                    title="Simulation pour essayer");



.. image:: 2020_covid_49_0.png


On compare avec les vraies données en gras.

.. code:: ipython3

    beta = 0.04
    mu = 0.03
    nu = 0.0001
    S0, I0, R0, D0 = france[120, :]
    sim = simulation_cumulee(beta, mu, nu, S0, I0, R0, D0, 120)
    plot_simulation(sim, dates[120], safe=False, two=True, true_data=france[120:240],
                    title="Simulation et vraies données (en gros)");



.. image:: 2020_covid_51_0.png


Pas simple de choisir des paramètres pour approximer la courbe.

.. code:: ipython3

    def error(data, simulation):
        err = (data[:, 1] - simulation[:, 1]) ** 2 + (data[:, 3] - simulation[:, 3]) ** 2
        total = (numpy.sum(err) / data.shape[0]) ** 0.5
        return total
    
    
    beta = 0.5
    mu = 1./14
    nu = 1./21
    S0 = 9990
    I0 = 10
    R0 = 0
    D0 = 0
    
    sim = simulation_cumulee(beta, mu, nu, S0, I0, R0, D0, 30)
    plot_simulation(sim);



.. image:: 2020_covid_53_0.png


On regarde si on arrive à retrouver les paramètres de la simulation.

Q3, Q4, Q5 : optimisation sur des données synthétiques
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. code:: ipython3

    from tqdm import tqdm  # pour avoir une barre de progression
    
    
    def optimisation(true_data, i_range=(0, 0.2), beta_range=(0, 0.5),
                     mu_range=(0., 0.2), nu_range=(0., 0.2),
                     max_iter=1000, error_fct=error):
        N = sum(true_data[0, :])
        rnd = numpy.random.rand(max_iter, 4)
        for i, (a, b) in enumerate([i_range, beta_range, mu_range, nu_range]):
            rnd[:, i] = rnd[:, i] * (b - a) + a
        err_min = None
        for it in tqdm(range(max_iter)):
            i, beta, mu, nu = rnd[it, :]
            dI0 = true_data[0, 0] * i
            D0 = true_data[0, 3]
            S0 = true_data[0, 0] - dI0
            I0 = true_data[0, 1] + dI0
            R0 = N - D0 - I0 - S0
            sim = simulation_cumulee(beta, mu, nu, S0, I0, R0, D0, days=true_data.shape[0] - 1)
            err = error_fct(true_data, sim)
            if err_min is None or err < err_min:
                err_min = err
                best = dict(beta=beta, mu=mu, nu=nu, I0=I0, i=i,
                            S0=S0, R0=R0, D0=D0, err=err, sim=sim)
        return best
    
    
    beta = 0.04
    mu = 0.07
    nu = 0.04
    S0 = 67e6
    I0 = 100000
    R0 = 10000
    D0 = 10000
    sim = simulation_cumulee(beta, mu, nu, S0, I0, R0, D0, 30)
    
    
    res = optimisation(sim, max_iter=2000, error_fct=error, i_range=(0., 0.001))
    sim_opt = res['sim']
    del res['sim']
    plot_simulation(sim_opt, dates[90], safe=False, two=True, true_data=sim,
                    title="beta=%1.3f mu=%1.3f nu=%1.3f err=%1.3g i=%1.3g" % (
                        res['beta'], res['mu'], res['nu'], res['err'], res['i']));


.. parsed-literal::
    100%|██████████| 2000/2000 [00:00<00:00, 7060.32it/s]



.. image:: 2020_covid_56_1.png


Ca ne marche pas trop mal.

Q6 : sur des données réelles
~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. code:: ipython3

    res = optimisation(france[180:], error_fct=error, i_range=(0, 0.0001))
    sim = res['sim']
    del res['sim']
    plot_simulation(sim, dates[180], safe=False, two=True, true_data=france[180:],
                    title="beta=%1.3f mu=%1.3f nu=%1.3f err=%1.3f i=%1.3f" % (
                        res['beta'], res['mu'], res['nu'], res['err'], res['i']));


.. parsed-literal::
    100%|██████████| 1000/1000 [00:03<00:00, 329.94it/s]



.. image:: 2020_covid_59_1.png


Ca ne marche pas super. On peut modifier l’erreur pour donner plus
d’importance à la courbe des morts.

.. code:: ipython3

    def error_norm(data, simulation):
        m1 = numpy.max(simulation[:, 1])
        m3 = numpy.max(simulation[:, 3])
        err = (data[:, 1] - simulation[:, 1]) ** 2 / m1 ** 2 + (data[:, 3] - simulation[:, 3]) ** 2 / m3 ** 2
        total = (numpy.sum(err) / data.shape[0]) ** 0.5
        return total
    
    res = optimisation(france[150:], error_fct=error_norm, i_range=(0, 0.001))
    sim = res['sim']
    del res['sim']
    plot_simulation(sim, dates[150], safe=False, two=True, true_data=france[150:],
                    title="beta=%1.3f mu=%1.3f nu=%1.3f err=%1.3f i=%1.3f" % (
                        res['beta'], res['mu'], res['nu'], res['err'], res['i']));


.. parsed-literal::
    100%|██████████| 1000/1000 [00:03<00:00, 314.34it/s]



.. image:: 2020_covid_61_1.png


L’erreur est répartie sur l’ensemble de la courbe mais cela ne tient pas
compte de la dynamique. Il faudrait prendre en compte les dérivées comme
dans l’aticle `A Modified SIR Model for the COVID-19 Contagion in
Italy <https://arxiv.org/pdf/2003.14391.pdf>`__. Ou partir sur une
approche plus directe comme celle proposée dans l’article `Estimating
and Simulating a SIRD Model of COVID-19 for Many Countries, States, and
Cities <https://web.stanford.edu/~chadj/sird-paper.pdf>`__.
:math:`\frac{dS}{dt} = - \beta \frac{S I}{N}` implique que
:math:`\Delta(S_t) = S_{t+1} - S_t = -\beta \frac{S_t I_t}{N}`. On peut
alors calculer une sorte de :math:`\beta_t` dépendant du temps :
:math:`\beta_t = -N\frac{S_{t+1} - S_t}{S_t I_t}`. Le problème, c’est
:math:`I_t` qu’on ne connaît pas puisque c’est le nombre de personnes
contaminantes à un instant *t*. On ne connaît pas plus :math:`S_t` mais
on peut le déduire si on connaît :math:`I_t`. Cela dit, le paramètre
:math:`\nu` est plus ou moins connu puisqu’il s’agit de l’inverse de la
durée moyenne d’incubation jusqu’au décès. De là, on peut retrouver
:math:`\beta`. On poursuit en passant au logarithme
:math:`\ln \beta_t = \ln N + \ln (S_t - S_{t+1}) - \ln S_t - \ln I_t`.
Le nombre de cas positifs découverts chaque jour correspond à
:math:`S_t - S_{t+1}`. On ajoute que
:math:`\forall t, N = S_t + I_t + R_t + D_t` et
:math:`D_{t+1} - D_t = \nu I_t`. Donc :

.. math:: \ln \beta_t = \ln N + \ln (S_t - S_{t+1}) - \ln (S_t - S_{t-1} + S_{t-1}) - \ln (D_{t+1} - D_t) + \ln \nu

Il y a deux inconnues, :math:`\ln \beta - \ln \nu` qu’on suppose
constant pendant la période et :math:`S_{t-1}`. Il suffit d’écrire ces
équations sur quelques jours puis de résoudre le système d’équations.

.. math::

   \begin{array}{l}
   \ln \beta - \ln \nu = \ln N + \ln (S_t - S_{t+1}) - \ln (S_t - S_{t-1} + S_{t-1}) - \ln (D_{t+1} - D_t) \\
   \ln \beta - \ln \nu = \ln N + \ln (S_{t+1} - S_{t+2}) - \ln (S_{t+1} - S_{t} + S_{t} - S_{t-1} + S_{t-1}) - \ln (D_{t+2} - D_{t+1}) 
   \end{array}

D’autres directions sont possibles comme `The Parameter Identification
Problem for SIR Epidemic Models: Identifying Unreported
Cases <https://www.math.u-bordeaux.fr/~pmagal100p/papers/MW-JMB-2018.pdf>`__
qui propose d’estimer les paramètres sur une plus grande période.

D’autres articles : `Introduction to inference: parameter
estimation <https://kingaa.github.io/clim-dis/parest/parest.html>`__,
`Data-based analysis, modelling and forecasting of the COVID-19
outbreak <https://europepmc.org/article/med/32231374>`__,
`Identification and estimation of the SEIRD epidemic model for
COVID-19 <https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7392128/>`__.

Quelques liens :


-  `CovidTracker <https://covidtracker.fr/covidtracker-france/>`__
-  `pyepidemics <https://github.com/collectif-codata/pyepidemics>`__
-  `COVID-19 pandemics
   control <https://freakonometrics.hypotheses.org/61300>`__
-  `modcov19 <https://modcov19.math.cnrs.fr/>`__
-  `aftercovid <http://www.xavierdupre.fr/app/aftercovid/helpsphinx/index.html>`__
-  `COVID-19 -
   France <https://dashboard.covid19.data.gouv.fr/vue-d-ensemble?location=FRA>`__
-  `Awesome COVID19
   resources <https://github.com/collectif-codata/awesome-covid19>`__