Modèle épidémiologique (SIRD)#

Etrangement, il n’était pas possible d’accéder aux données de l’épidémie le jour du TD car le site data.gouv.fr était inacessible suite à une panne chez OVH.

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

Comme les vraies données sont inaccessibles, on s’intéresse à des données simulées.

Simulation#

[3]:

import numpy


def iteration(S, I, R, D, beta, mu, nu, N):
    dR = I * mu
    dD = I * nu
    dS = -S * I * beta / N
    dI = -dS - dR - dD
    return S + dS, I + dI, R + dR, D + dD


iteration(10, 1, 0, 0, 0.1, 0.1, 0.1, 11)

[3]:

(9.909090909090908, 0.8909090909090909, 0.1, 0.1)

[4]:

def simulation(n_iter, S, I, R, D, beta, mu, nu):
    sim = numpy.zeros((n_iter + 1, 4))
    sim[0, :] = numpy.array([S, I, R, D]).ravel()
    N = S + I + R + D
    for i in range(n_iter):
        S, I, R, D = iteration(S, I, R, D, beta, mu, nu, N)
        sim[i + 1, :] = numpy.array([S, I, R, D]).ravel()
    return sim


simulation(2, 10, 1, 0, 0, 0.1, 0.1, 0.1)

[4]:

array([[10.        ,  1.        ,  0.        ,  0.        ],
       [ 9.90909091,  0.89090909,  0.1       ,  0.1       ],
       [ 9.82883546,  0.79298272,  0.18909091,  0.18909091]])

Visualisation#

[5]:

%matplotlib inline

[7]:

import matplotlib.pyplot as plt
import pandas

fig, ax = plt.subplots(1, 1, figsize=(10, 2))
sim = simulation(20, 10, 0.5, 0, 0, 1.2, 0.2, 0.1)
df = pandas.DataFrame(sim, columns=["S", "I", "R", "D"])
df.plot(ax=ax);

../../_images/practice_algo-base_exercice_simulation_covid_simple_7_0.png

Estimation des paramètres#

Une fois qu’on a des données, on peut essayer de retrouver la paramètres de la simulation. Mais auparavant, on simule des données un peu plus bruitées en changeant les paramètres du modèle tous les jours.

[8]:

def simulation_bruitee(n_iter, S, I, R, D, beta, mu, nu):
    sim = numpy.zeros((n_iter + 1, 4))
    sim[0, :] = numpy.array([S, I, R, D]).ravel()
    N = S + I + R + D
    for i in range(n_iter):
        b = numpy.random.randn(1) * beta / 5 + beta
        m = numpy.random.randn(1) * mu / 5 + mu
        n = numpy.random.randn(1) * nu / 5 + nu
        S, I, R, D = iteration(S, I, R, D, b, m, n, N)
        sim[i + 1, :] = numpy.array([S, I, R, D]).ravel()
    return sim


simulation_bruitee(2, 10, 1, 0, 0, 0.1, 0.1, 0.1)

[8]:

array([[10.        ,  1.        ,  0.        ,  0.        ],
       [ 9.90629495,  0.90818313,  0.0769856 ,  0.10853633],
       [ 9.84322669,  0.77113049,  0.17022418,  0.21541865]])

[9]:

fig, ax = plt.subplots(1, 1, figsize=(10, 2))
sim_bruit = simulation_bruitee(20, 10, 0.5, 0, 0, 1.2, 0.2, 0.1)
df = pandas.DataFrame(sim_bruit, columns=["S", "I", "R", "D"])
df.plot(ax=ax);

../../_images/practice_algo-base_exercice_simulation_covid_simple_10_0.png

Première idée#

On calcule les paramètres chaque jour à partir de la définition du modèle.

[10]:

def estimation_coefficient(sim):
    N = sim.sum(axis=1)
    diff = sim[1:] - sim[:-1]
    mu = diff[:, 2] / sim[:-1, 1]
    nu = diff[:, 3] / sim[:-1, 1]
    beta = -diff[:, 0] / (sim[:-1, 1] * sim[:-1, 0]) * N[:-1]
    return dict(beta=beta, mu=mu, nu=nu)


df = pandas.DataFrame(estimation_coefficient(sim_bruit))
df.head(n=1)

[10]:

	beta	mu	nu
0	0.645167	0.153483	0.113224

[11]:

df.plot(figsize=(10, 2));

../../_images/practice_algo-base_exercice_simulation_covid_simple_13_0.png

On pourrait imaginer que les paramètres cherchés soient égaux à la moyenne des valeurs obtenues pour chaque jour. Encore faudrait-il le prouver ou trouver un contre exemple.

Seconde idée : problèmes d’optimisation#

On cherche à transformer le problème d’estimation en un problème d’optimisation de telle sorte que les paramètres en sont la solution. L’idée de simuler en fonction d’un jeu de paramètres puis de calculer une distance entre cette simulation et les données observées.

Le problème de cette solution est qu’il faut essayer plein de valeur, au hasard pour faire simple. Donc on se sert de la fonction précédente comme point de départ et on tire des valeurs tout autour.

[12]:

def distance_sim(obs, sim):
    return ((obs - sim) ** 2).sum() / obs.size


distance_sim(sim, sim_bruit)

[12]:

0.1379327480825296

[13]:

def optimisation(obs):
    est = estimation_coefficient(obs)
    b0 = est["beta"].mean()
    m0 = est["mu"].mean()
    n0 = est["nu"].mean()
    n_iter = obs.shape[0] - 1
    S, I, R, D = obs[0]
    beta, mu, nu = b0, m0, n0
    distance = distance_sim(obs, simulation(n_iter, S, I, R, D, b0, m0, n0))
    distance0 = distance
    for i in range(0, 100):
        b = numpy.random.randn(1) * b0 / 5 + b0
        m = numpy.random.randn(1) * m0 / 5 + m0
        n = numpy.random.randn(1) * n0 / 5 + n0
        sim = simulation(n_iter, S, I, R, D, b, m, n)
        d = distance_sim(obs, sim)
        if distance is None or d < distance:
            distance = d
            beta, mu, nu = b, m, n
    return dict(
        b0=b0,
        mu0=m0,
        nu0=n0,
        beta=beta,
        mu=mu,
        nu=nu,
        distance=distance,
        distance0=distance0,
    )


optimisation(sim_bruit)

[13]:

{'b0': 1.175492030437819,
 'mu0': 0.21093122422461147,
 'nu0': 0.11141742761811038,
 'beta': array([1.12240477]),
 'mu': array([0.1995322]),
 'nu': array([0.09668045]),
 'distance': 0.06764482535407212,
 'distance0': 0.0880922396685696}

distance0 la distance entre la simulation obtenue avec les paramètes de la première méthode, distance celle de la seconde et elle est plus petite. Cela répond à la question qu’on se posait ci-dessus.

[ ]:

Notebook on github