[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

Essayez le langage de programmation stochastique Stan

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

«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.

Les données des personnes infectées symptomatiques datent de 2/20 et le nombre d'infections a encore augmenté lorsque le 7e décès a été confirmé (3/7), mais même en considérant que, COVID Il est peu probable que les chances de mourir dépassent 1%.

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);
}

Stan peut gérer les tableaux, donc je peux écrire plus clairement, mais je n'aimais pas laisser tomber la classe dans les indices du tableau, donc je les énumère honnêtement.

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.

Si vous prenez simplement floor, le dégradé disparaîtra et il ne convergera pas bien, vous devez donc bien le lisser en utilisant ʻinv_logit` ou quelque chose.

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

https://www.jiji.com/jc/article?k=2020030501355&g=soc ↩︎
https://www.niid.go.jp/niid/ja/diseases/ka/corona-virus/2019-ncov/2484-idsc/9422-covid-dp-2.html ↩︎ ↩︎
https://www.yomiuri.co.jp/national/20200307-OYT1T50228/ ↩︎

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

Objectif

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

Préparation

Les données

Environnement d'exécution

Pour le moment de Jab

Tenez compte de l'âge

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

Vérifiez d'abord avec NumPy

Vérifions avec Stan

Essayez de résoudre avec la distribution binomiale

Essayez de résoudre avec la distribution bêta

Résumé

appendice