La nouvelle Corona est-elle vraiment une menace? Validé avec Stan (était)

[Avis] Sequel demande le taux de létalité du COVID-19.

Objectif

--Vérifiez si la létalité du COVID-19 est vraiment élevée

Qu'est-ce qu'un langage de programmation probabiliste?

C'est un langage de programmation qui résout les relations statistiques (modèles) d'une manière agréable. En tant que mise en œuvre, il semble effectuer un échantillonnage par la méthode de Monte Carlo en chaîne de Markov.

Il ne remplace pas les langages de programmation procédurale car ils ont des rôles et des méthodes complètement différents.

Stan est l'un des langages de programmation probabilistes.

Préparation

Les données

Utilisez le statut d'infection du Diamond Princess. En effet, il n'y a pas d'espace dans le monde où le statut de l'infection a été étudié plus en détail, et on s'attend à ce que la vérification la plus précise soit possible.

article Les données Date et l'heure La source
Nombre de personnes infectées 696 personnes 3/5 Affaires courantes dot com[1]
Nombre de personnes infectées(Symptôme) 301 personnes 2/20 Institut national des maladies infectieuses[2]
Nombre de personnes infectées(Asymptomatique) 318 personnes 2/20 Institut national des maladies infectieuses[2:1]
Nombre de décès 7 personnes 3/7 Yomiuri Shimbun[3]

Nous voulons également comparer avec la grippe par âge, nous allons donc traiter les statistiques du CDC [^ 4] afin de pouvoir les comparer avec les données de l'Institut national des maladies infectieuses [^ 2]. La conversion de classe a été faite simplement par son rapport de largeur. Étant donné que les statistiques du CDC ne contiennent que des chiffres sur le nombre d'infections grippales symptomatiques, nous utiliserons également les données sur les infections symptomatiques pour le COVID-19.

classe CODIV-19 Nombre d'infections symptomatiques Létalité grippale(Estimation)
0-9 0 0.0050%
10-19 2 0.0063%
20-29 25 0.0206%
30-39 27 0.0206%
40-49 19 0.0206%
50-59 28 0.0614%
60-69 76 0.4465%
70-79 95 0.8315%
80-89 27 0.8315%
90-99 2 0.8315%
Total 301 0.0962%

Environnement d'exécution

Installez Stan. Stan peut être exécuté uniquement sur la ligne de commande, mais il est beaucoup plus facile d'utiliser un wrapper.

Cette fois, nous utiliserons PyStan, une interface Python qui peut être facilement introduite avec pip. Vous avez besoin de numpy et cython, ainsi que de scipy et matplotlib pour afficher le graphique, donc incluez-les également.

Lors de l'utilisation d'Anaconda

$ conda create -n dp-corona-stan python=3.7 numpy cython scipy matplotlib pystan
$ conda activate dp-corona-stan

Pour le moment de Jab

La létalité est estimée par Stan à partir du nombre de personnes infectées (696 personnes) et du nombre de décès (7 personnes).

import pystan

model = pystan.StanModel(model_code='''
data {
    int N; //Nombre de personnes infectées
    int D; //Nombre de décès
}
parameters {
    real<lower=0, upper=1> p;
}
model {
    D ~ binomial(N, p);
}
''')

data = {
    'N': 696,
    'D':   7
}

fit = model.sampling(data=data, chains=4)
print(fit)

La chaîne passée à StanModel est le code Stan.

Cet article ne donne pas trop de détails sur la façon d'écrire Stan, mais en un mot, cela signifie entrer les données spécifiées dans data et optimiser les variables spécifiées dans parameters pour satisfaire model. C'est un code.

L'événement du nombre de personnes infectées décédées suivra statistiquement une distribution binomiale. Autrement dit, la description de «modèle»,

    D ~ binomial(N, p);

Est un modèle de la situation où des personnes «D» sont décédées lorsque chaque personne infectée «N» est décédée avec une probabilité «p».

Vous pouvez exécuter ce code pour estimer le taux de létalité «p» pour toutes les personnes infectées.

       mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
p      0.01  9.8e-5 4.1e-3 5.0e-3 8.4e-3   0.01   0.01   0.02   1724    1.0
lp__ -44.22    0.02   0.73  -46.3  -44.4 -43.95 -43.76  -43.7   1531    1.0

corona_death_estimation_stan_simple.png

«p» est estimé à environ 1%. Le public dit qu'il est de 2%, mais à partir de ces données, 2% sont sur la ligne de coup de l'intervalle de confiance à 95% («5,0e-3» à «0,02»), et la possibilité est rejetée. Je ne peux pas deviner, mais ce ne sera pas si cher.

Cependant, la moitié des personnes infectées sont asymptomatiques, donc si nous repositionnons la létalité et les critères pour les personnes infectées symptomatiques, ce sera certainement autour de 2%. (À propos, il semble que 1/3 de la grippe soit asymptomatique [^ 5].)

Tenez compte de l'âge

Cependant, plus de la moitié de l'équipage et des passagers du Diamond Princess avaient plus de 60 ans. Ceci est différent de la répartition de la population dans le monde, il est donc nécessaire d'en tenir compte pour obtenir un taux de mortalité réaliste.

Par conséquent, si le taux de létalité du COVID-19 est équivalent à celui de la grippe, nous estimons le nombre de décès qui se produiront en cas d'infection à l'échelle de la princesse du diamant et le comparons au nombre réel de décès (7 personnes). Je vais.

Préparation: représenter les données en Python

Enregistrez les données préparées dans Data dans un fichier py.

«S_xx» est le nombre estimé d'infections symptomatiques par groupe d'âge, et «p_flu_xx» est le taux de mortalité estimé pour le nombre d'infections grippales symptomatiques par groupe d'âge.

# data.py
data = {
    'S_0x':  0,
    'S_1x':  2,
    'S_2x': 25,
    'S_3x': 27,
    'S_4x': 19,
    'S_5x': 28,
    'S_6x': 76,
    'S_7x': 95,
    'S_8x': 27,
    'S_9x':  2,
    'p_flu_0x': 0.000050,
    'p_flu_1x': 0.000063,
    'p_flu_2x': 0.000206,
    'p_flu_3x': 0.000206,
    'p_flu_4x': 0.000206,
    'p_flu_5x': 0.000614,
    'p_flu_6x': 0.004465,
    'p_flu_6x': 0.008315,
    'p_flu_7x': 0.008315,
    'p_flu_8x': 0.008315,
    'p_flu_9x': 0.000962
}

Vérifiez d'abord avec NumPy

NumPy a la capacité de générer des nombres aléatoires qui suivent une distribution binomiale. Si vous utilisez cette fonction, vous pouvez facilement effectuer des simulations.

import numpy as np
from data import data

N = 10000
MIN_DEATH = 7

sample = (np.random.binomial(data['S_0x'], data['p_flu_0x'], N) +
          np.random.binomial(data['S_1x'], data['p_flu_1x'], N) +
          np.random.binomial(data['S_2x'], data['p_flu_2x'], N) +
          np.random.binomial(data['S_3x'], data['p_flu_3x'], N) +
          np.random.binomial(data['S_4x'], data['p_flu_4x'], N) +
          np.random.binomial(data['S_5x'], data['p_flu_5x'], N) +
          np.random.binomial(data['S_6x'], data['p_flu_6x'], N) +
          np.random.binomial(data['S_7x'], data['p_flu_7x'], N) +
          np.random.binomial(data['S_8x'], data['p_flu_8x'], N) +
          np.random.binomial(data['S_9x'], data['p_flu_9x'], N))
probability = float(sum(sample >= MIN_DEATH)) / N

print('Average # of deaths: %.2f' % (sample.mean()))
print('# of deaths >= %d: %.2f%%' % (MIN_DEATH, probability * 100))

Je n'entrerai pas dans cela aussi, mais je vais l'expliquer brièvement.

np.random.binomial (n, p, N) répète la tâche d'essayer les événements qui se produisent avec une probabilité p n fois et enregistre le nombre d'occurrences N fois. Par exemple, np.random.binomial (2, 0.5, 3) retournera ʻarray ([2, 0, 1]) `.

Un tableau de NumPy est un ajout élément par élément lorsqu'il est ajouté normalement, ce qui donne le nombre d'occurrences dans toutes les classes.

De plus, si vous utilisez l'opérateur de comparaison sur un tableau NumPy, il renverra le tableau créé par le jugement d'authenticité pour chaque élément. Par exemple, np.array ([2, 0, 1])> 1 renverra ʻarray ([True, False, False])`.

La fonction «sum» de Python considère «True» comme «1» et «False» comme «0», donc après tout, «float (sum (sample> = MIN_DEATH)) / N» est utilisé pour essayer «N». , Vous pouvez calculer le pourcentage de décès dépassant «MIN_DEATH».

Le résultat est

Average # of deaths: 1.67
# of deaths >= 7: 0.21%

C'est devenu comme.

Si le taux de mortalité du COVID-19 est similaire à celui de la grippe, il semble tout à fait improbable que jusqu'à 7 décès surviennent pour 301 personnes infectées symptomatiques.

Vérifions avec Stan

Essayez de résoudre avec la distribution binomiale

Vous devriez être en mesure d'écrire quelque chose de similaire au code Stan ci-dessus, la seule différence étant que la cible d'optimisation passe de la probabilité aux décès et qu'elle prend en compte le groupe d'âge.

Si le nombre estimé de décès par âge est «d_xx»,

data {
    int S_0x;
    int S_1x;
    int S_2x;
    int S_3x;
    int S_4x;
    int S_5x;
    int S_6x;
    int S_7x;
    int S_8x;
    int S_9x;
    real p_flu_0x;
    real p_flu_1x;
    real p_flu_2x;
    real p_flu_3x;
    real p_flu_4x;
    real p_flu_5x;
    real p_flu_6x;
    real p_flu_7x;
    real p_flu_8x;
    real p_flu_9x;
}
parameters {
    // upper = S_xx + 1 so that S_xx can be 0
    int<lower=0, upper=S_0x+1> d_0x;
    int<lower=0, upper=S_1x+1> d_1x;
    int<lower=0, upper=S_2x+1> d_2x;
    int<lower=0, upper=S_3x+1> d_3x;
    int<lower=0, upper=S_4x+1> d_4x;
    int<lower=0, upper=S_5x+1> d_5x;
    int<lower=0, upper=S_6x+1> d_6x;
    int<lower=0, upper=S_7x+1> d_7x;
    int<lower=0, upper=S_8x+1> d_8x;
    int<lower=0, upper=S_9x+1> d_9x;
}
transformed parameters {
    int d;
    d = d_0x + d_1x + d_2x + d_3x + d_4x + d_5x + d_6x + d_7x + d_8x + d_9x;
}
model {
    d_0x ~ binomial(S_0x, p_flu_0x);
    d_1x ~ binomial(S_1x, p_flu_1x);
    d_2x ~ binomial(S_2x, p_flu_2x);
    d_3x ~ binomial(S_3x, p_flu_3x);
    d_4x ~ binomial(S_4x, p_flu_4x);
    d_5x ~ binomial(S_5x, p_flu_5x);
    d_6x ~ binomial(S_6x, p_flu_6x);
    d_7x ~ binomial(S_7x, p_flu_7x);
    d_8x ~ binomial(S_8x, p_flu_8x);
    d_9x ~ binomial(S_9x, p_flu_9x);
}

Quand tu cours

ValueError: Failed to parse Stan model 'anon_model_fecb1e77228fe372ef4eb9bc4bcc8086'. Error message:
SYNTAX ERROR, MESSAGE(S) FROM PARSER:
Parameters or transformed parameters cannot be integer or integer array;  found int variable declaration, name=d_0x
 error in 'unknown file name' at line 26, column 35
  -------------------------------------------------
    24: parameters {
    25:     // upper = S_xx + 1 so that S_xx can be 0
    26:     int<lower=0, upper=S_0x+1> d_0x;
                                          ^
    27:     int<lower=0, upper=S_1x+1> d_1x;
  -------------------------------------------------

J'étais en colère parce que je ne pouvais pas utiliser un entier. .. ..

À la page 20 du groupe de lecture de modélisation statistique bayésienne Ch.9 [^ 6] avec Stan et R, la faiblesse de Stan était que ʻintne pouvait pas être utilisé pour lesparamètres`. Je ne peux pas vraiment m'en servir.

Ω \ ζ ゜) Chaîne ♪

Puisque binomial demande ʻint sur le côté gauche, l'astuce de déclarer d_xxs avec real au lieu de ʻint ne peut pas être utilisée.

Ω \ ζ ゜) Chaîne ♪

Essayez de résoudre avec la distribution bêta

La distribution binomiale ne fonctionnait pas, j'ai donc essayé d'utiliser la distribution bêta.

En parlant de la préparation que Masakari volera, la distribution bêta est l'opposé de la distribution binomiale, qui est la distribution de probabilité de la probabilité d'occurrence «p» dans le cas qui suit la distribution binomiale. Cliquez ici pour plus de détails.

    d_xx ~ binomial(S_xx, p_flu_xx);

Ce qui a été écrit

    p_flu_xx ~ beta(d_xx + 1, S_xx - d_xx + 1);

Réécrivez et changez le type de «d_xx» en «réel».

import pystan
from data import data

model = pystan.StanModel(model_code="""
data {
    int S_0x;
    int S_1x;
    int S_2x;
    int S_3x;
    int S_4x;
    int S_5x;
    int S_6x;
    int S_7x;
    int S_8x;
    int S_9x;
    real p_flu_0x;
    real p_flu_1x;
    real p_flu_2x;
    real p_flu_3x;
    real p_flu_4x;
    real p_flu_5x;
    real p_flu_6x;
    real p_flu_7x;
    real p_flu_8x;
    real p_flu_9x;
}
parameters {
    // upper = S_xx + 1 so that S_xx can be 0
    real<lower=0, upper=S_0x+1> d_0x;
    real<lower=0, upper=S_1x+1> d_1x;
    real<lower=0, upper=S_2x+1> d_2x;
    real<lower=0, upper=S_3x+1> d_3x;
    real<lower=0, upper=S_4x+1> d_4x;
    real<lower=0, upper=S_5x+1> d_5x;
    real<lower=0, upper=S_6x+1> d_6x;
    real<lower=0, upper=S_7x+1> d_7x;
    real<lower=0, upper=S_8x+1> d_8x;
    real<lower=0, upper=S_9x+1> d_9x;
}
transformed parameters {
    real d;
    d = d_0x + d_1x + d_2x + d_3x + d_4x + d_5x + d_6x + d_7x + d_8x + d_9x;
}
model {
    p_flu_0x ~ beta(d_0x + 1, S_0x - d_0x + 1);
    p_flu_1x ~ beta(d_1x + 1, S_1x - d_1x + 1);
    p_flu_2x ~ beta(d_2x + 1, S_2x - d_2x + 1);
    p_flu_3x ~ beta(d_3x + 1, S_3x - d_3x + 1);
    p_flu_4x ~ beta(d_4x + 1, S_4x - d_4x + 1);
    p_flu_5x ~ beta(d_5x + 1, S_5x - d_5x + 1);
    p_flu_6x ~ beta(d_6x + 1, S_6x - d_6x + 1);
    p_flu_7x ~ beta(d_7x + 1, S_7x - d_7x + 1);
    p_flu_8x ~ beta(d_8x + 1, S_8x - d_8x + 1);
    p_flu_9x ~ beta(d_9x + 1, S_9x - d_9x + 1);
}
""")

fit = model.sampling(data=data, chains=4)
print(fit)

Voici le résultat de l'exécution:

       mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
d_0x    0.1  1.4e-3   0.09 2.3e-3   0.03   0.07   0.14   0.35   4361    1.0
d_1x   0.12  1.7e-3   0.12 2.5e-3   0.03   0.08   0.16   0.43   4474    1.0
d_2x   0.19  2.6e-3   0.18 5.0e-3   0.06   0.14   0.27   0.67   4976    1.0
d_3x    0.2  3.1e-3   0.19 4.4e-3   0.06   0.14   0.27   0.71   3641    1.0
d_4x   0.19  2.5e-3   0.19 3.5e-3   0.05   0.13   0.26   0.69   5315    1.0
d_5x   0.25  3.2e-3   0.23 8.0e-3   0.08   0.18   0.34   0.85   5100    1.0
d_6x   0.93    0.01   0.73   0.04   0.36   0.77   1.33    2.8   4387    1.0
d_7x   1.06    0.01    0.8   0.04   0.43    0.9   1.51   3.01   4028    1.0
d_8x   0.54  7.0e-3   0.48   0.01   0.18   0.42   0.76   1.82   4729    1.0
d_9x   0.17  2.4e-3   0.16 4.3e-3   0.05   0.12   0.24   0.58   4580    1.0
d      3.73    0.02   1.25   1.66   2.84   3.59   4.51   6.57   4367    1.0
lp__  -1.64    0.08   2.58  -7.78  -3.13  -1.28   0.24   2.37   1048    1.0

«d» est le nombre de décès, mais «50%» est «3,59», ce qui est nettement supérieur au «1,67» obtenu par NumPy.

Ω \ ζ ゜) Chaîne ♪

Je me demande pourquoi.

C'est probablement un problème car «d» est défini comme «réel». Il semble que d_xx, qui à l'origine ne peut prendre que des valeurs entières, peut prendre des valeurs fractionnaires, ce qui s'écarte du problème réel.

En fait, si vous essayez de calculer en convertissant d_xx en un entier avec étage

       mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
...
d      1.82    0.03   1.29   0.01   1.01   1.69    2.5   4.93   1710    1.0

C'est devenu une valeur comme ça.

Cependant, si vous souhaitez faire cela, vous pouvez utiliser NumPy. Il est difficile de dire que c'est plus correct pour la modélisation stochastique.

Existe-t-il un bon moyen de résoudre ce problème avec Stan?

Résumé

En partie, nous avons pu utiliser Stan pour estimer le risque du nouveau virus corona.

D'après la situation du Diamond Princess, il semble tout à fait certain que le COVID-19 est plus vicieux en termes de létalité que la grippe traditionnelle. Une autre menace est que si de nombreuses personnes ont déjà des anticorps contre la grippe, peu d'entre elles ont des anticorps COVID-19. Vous devez être vigilant pour le bien de la société dans son ensemble, plutôt que pour l'individu.

Cependant, on estime qu'environ 1,67 personne mourra si elle est infectée par la grippe à cette échelle, donc si vous avez tellement peur du COVID-19, vous devriez être plus prudent avec la grippe ordinaire. De mon point de vue personnel, je voudrais faire une déclaration.

Sequel conteste également la tâche de trouver la létalité du COVID-19. S'il vous plaît, jetez un oeil.

appendice

Code source: https://github.com/akeyhero/dp-corona-stan


  1. https://www.jiji.com/jc/article?k=2020030501355&g=soc ↩︎

  2. https://www.niid.go.jp/niid/ja/diseases/ka/corona-virus/2019-ncov/2484-idsc/9422-covid-dp-2.html ↩︎ ↩︎

  3. https://www.yomiuri.co.jp/national/20200307-OYT1T50228/ ↩︎

Recommended Posts

La nouvelle Corona est-elle vraiment une menace? Validé avec Stan (était)
[Suite] La nouvelle Corona est-elle vraiment une menace? Létalité calculée avec Stan
[New Corona] Le prochain pic est-il en décembre? J'ai essayé l'analyse des tendances avec Python!
Créer un nouveau csv avec des pandas basé sur le csv local
L'image est Namekuji
J'ai essayé de prédire le nombre de personnes infectées au niveau national de la nouvelle corona avec un modèle mathématique
Essayez de jouer avec un simulateur de baseball # 1 La carie d'alimentation est-elle efficace?
Erreur avec pip: un problème est survenu lors de la confirmation du certificat SSL
[Python] Qu'est-ce qu'une instruction with?
L'espace est dangereux avec PyEphem
Création d'une nouvelle application corona à Kyoto avec le framework Web de Python Dash
Le panneau Web LXC qui peut faire fonctionner LXC avec un navigateur était merveilleux
J'ai essayé de créer facilement une image 3D de haute précision avec une seule photo [-1]. (La zone cachée est-elle vraiment visible?)