[PYTHON] Ist die neue Corona wirklich eine Bedrohung? Validiert mit Stan (war)

[Hinweis] Fortsetzung fragt nach der Letalitätsrate von COVID-19.

Zweck

Was ist eine probabilistische Programmiersprache?

Es ist eine Programmiersprache, die statistische Beziehungen (Modelle) auf nette Weise löst. Als Implementierung scheint es eine Stichprobe nach der Markov-Ketten-Monte-Carlo-Methode durchzuführen.

Es ist kein Ersatz für prozedurale Programmiersprachen, da sie völlig unterschiedliche Rollen und Methoden haben.

Stan ist eine der probabilistischen Programmiersprachen.

Vorbereitung

Daten

Verwenden Sie den Infektionsstatus der Diamond Princess. Dies liegt daran, dass es auf der Welt keinen Raum gibt, in dem der Infektionsstatus genauer untersucht wurde, und es wird erwartet, dass die genaueste Überprüfung möglich ist.

Artikel Daten Datum (und Uhrzeit Quelle
Anzahl der infizierten Personen 696 Personen 3/5 Aktuelle Angelegenheiten dot com[1]
Anzahl der infizierten Personen(Symptom) 301 Personen 2/20 Nationales Institut für Infektionskrankheiten[2]
Anzahl der infizierten Personen(Asymptomatisch) 318 Personen 2/20 Nationales Institut für Infektionskrankheiten[2:1]
Anzahl der Todesfälle 7 Personen 3/7 Yomiuri Shimbun[3]

Ich möchte auch nach Alter mit Influenza vergleichen, daher werde ich die CDC-Statistiken [^ 4] verarbeiten, damit ich sie mit den Daten des Nationalen Instituts für Infektionskrankheiten [^ 2] vergleichen kann. Die Klassenkonvertierung erfolgte einfach anhand des Breitenverhältnisses. Da die CDC-Statistiken nur Zahlen zur Anzahl der symptomatischen Influenza-Infektionen enthalten, werden wir auch die Daten zu symptomatischen Infektionen für COVID-19 verwenden.

Klasse CODIV-19 Anzahl der symptomatischen Infektionen Influenza-Letalität(Schätzen)
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%

Ausführungsumgebung

Installieren Sie Stan. Stan kann allein über die Befehlszeile ausgeführt werden, es ist jedoch viel einfacher, einen Wrapper zu verwenden.

Dieses Mal werden wir PyStan verwenden, eine Python-Schnittstelle, die einfach mit pip eingeführt werden kann. Sie benötigen "numpy" und "cython" sowie "scipy" und "matplotlib", um das Diagramm anzuzeigen. Fügen Sie sie also ebenfalls hinzu.

Bei Verwendung von Anaconda

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

Vorerst von Jab

Die Letalität wird von Stan anhand der Anzahl der Infizierten (696 Personen) und der Anzahl der Todesfälle (7 Personen) geschätzt.

import pystan

model = pystan.StanModel(model_code='''
data {
    int N; //Anzahl der infizierten Personen
    int D; //Anzahl der Todesfälle
}
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)

Die Zeichenfolge, die Sie an "StanModel" übergeben, ist der Stan-Code.

Dieser Artikel geht nicht zu detailliert auf das Schreiben von Stan ein, aber kurz gesagt bedeutet dies, dass Sie die in "Daten" angegebenen Daten eingeben und die in "Parameter" angegebenen Variablen optimieren, um "Modell" zu erfüllen. Es ist ein Code.

Das Ereignis, wie viele der infizierten Menschen starben, folgt statistisch einer Binomialverteilung. Das heißt, die Beschreibung von "Modell",

    D ~ binomial(N, p);

Ist ein Modell der Situation, in der "D" -Personen starben, als jede infizierte Person "N" mit einer "p" -Wahrscheinlichkeit starb.

Sie können diesen Code ausführen, um die Letalitätsrate "p" für alle infizierten Personen zu schätzen.

       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" wird auf etwa 1% geschätzt. Es wird von der Öffentlichkeit gesagt, dass es 2% sind, aber von diesen Daten sind 2% auf der Trefferlinie des 95% -Konfidenzintervalls ("5,0e-3" bis "0,02"), und die Möglichkeit wird abgelehnt. Ich kann es nicht erraten, aber es wird nicht so teuer sein.

Die Hälfte der Infizierten ist jedoch asymptomatisch. Wenn wir also die Letalität und die Kriterien für symptomatische Infizierte neu positionieren, werden es sicherlich etwa 2% sein. (Übrigens scheint 1/3 der Influenza asymptomatisch zu sein [^ 5].)

Betrachten Sie das Alter

Mehr als die Hälfte der Besatzung und Passagiere von Diamond Princess waren jedoch über 60 Jahre alt. Dies unterscheidet sich von der Bevölkerungsverteilung in der Welt, daher muss dies berücksichtigt werden, um eine realistische Sterblichkeitsrate zu erhalten.

Wenn die Letalitätsrate von COVID-19 der von Influenza entspricht, schätzen wir daher, wie viele Todesfälle bei einer Infektion mit Diamantprinzessin auftreten werden, und vergleichen sie mit der tatsächlichen Anzahl der Todesfälle (7 Personen). Ich werde.

Vorbereitung: Stellen Sie die Daten in Python dar

Speichern Sie die in Daten vorbereiteten Daten in einer py-Datei.

"S_xx" ist die geschätzte Anzahl symptomatischer Infektionen nach Altersgruppen, und "p_flu_xx" ist die geschätzte Sterblichkeitsrate für die Anzahl symptomatischer Influenza-Infektionen nach Altersgruppen.

# 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
}

Überprüfen Sie dies zuerst mit NumPy

NumPy kann Zufallszahlen generieren, die einer Binomialverteilung folgen. Wenn Sie diese Funktion verwenden, können Sie problemlos Simulationen durchführen.

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

Ich werde auch nicht darauf eingehen, aber ich werde es kurz erklären.

np.random.binomial (n, p, N) wiederholt die Aufgabe, Ereignisse zu versuchen, die mit der Wahrscheinlichkeit p n mal auftreten, und die Anzahl der Vorkommen N mal aufzuzeichnen. Beispiel: "np.random.binomial (2, 0.5, 3)" gibt "array ([2, 0, 1])" zurück.

Ein Array von NumPy ist eine Element-für-Element-Addition, wenn es normal hinzugefügt wird. Dies gibt die Anzahl der Vorkommen in allen Klassen an.

Wenn Sie den Vergleichsoperator für ein NumPy-Array verwenden, wird das durch Authentizitätsbeurteilung für jedes Element erstellte Array zurückgegeben. Zum Beispiel wird "np.array ([2, 0, 1])> 1" ein Array ([True, False, False]) "zurückgeben.

Pythons "sum" -Funktion betrachtet "True" als "1" und "False" als "0". Schließlich wird "float (sum (sample> = MIN_DEATH)) / N" verwendet, um "N" zu versuchen. , Sie können den Prozentsatz der Todesfälle berechnen, die "MIN_DEATH" überschreiten.

Ergebnis ist

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

Es wurde wie.

Wenn die Sterblichkeitsrate von COVID-19 der von Influenza ähnlich ist, ist es ziemlich unwahrscheinlich, dass bei 301 symptomatisch infizierten Personen bis zu 7 Todesfälle auftreten.

Lassen Sie uns mit Stan überprüfen

Versuchen Sie, mit Binomialverteilung zu lösen

Sie sollten in der Lage sein, etwas Ähnliches wie den obigen Stan-Code zu schreiben, mit dem einzigen Unterschied, dass sich das Optimierungsziel von Wahrscheinlichkeit zu Tod ändert und die Altersgruppe berücksichtigt.

Wenn die geschätzte Anzahl der Todesfälle nach Alter "d_xx" ist,

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

Wenn du rennst

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

Ich war wütend, weil ich keine ganze Zahl verwenden konnte. .. ..

Auf Seite 20 der Bayesian Statistical Modeling Reading Group Ch.9 [^ 6] mit Stan und R war Stans Schwäche, dass "int" nicht für "Parameter" verwendet werden konnte. Ich kann es nicht wirklich benutzen.

Ω \ ζ ゜) Kette ♪

Da "binomial" auf der linken Seite nach "int" fragt, kann der Trick, "d_xx" mit "real" anstelle von "int" zu deklarieren, nicht verwendet werden.

Ω \ ζ ゜) Kette ♪

Versuchen Sie mit Beta-Distribution zu lösen

Die Binomialverteilung hat nicht funktioniert, deshalb habe ich versucht, die Beta-Verteilung zu verwenden.

In Bezug auf die Bereitschaft, dass Masakari fliegen wird, ist die Beta-Verteilung das Gegenteil der Binomialverteilung, die die Wahrscheinlichkeitsverteilung der Auftrittswahrscheinlichkeit "p" für den Fall ist, der auf die Binomialverteilung folgt. Klicken Sie hier für Details.

    d_xx ~ binomial(S_xx, p_flu_xx);

Was wurde geschrieben

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

Schreiben Sie um und ändern Sie den Typ von "d_xx" in "real".

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)

Hier sind die Ergebnisse des Laufens:

       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" ist die Anzahl der Todesfälle, aber "50%" ist "3,59", was deutlich höher ist als die von NumPy erhaltenen "1,67".

Ω \ ζ ゜) Kette ♪

Ich wundere mich warum.

Dies ist wahrscheinlich ein Problem, da "d" als "real" definiert ist. Es scheint, dass d_xx, das ursprünglich nur ganzzahlige Werte annehmen kann, Bruchwerte annehmen kann, die vom eigentlichen Problem abweichen.

In der Tat, wenn Sie versuchen zu berechnen, indem Sie "d_xx" in eine Ganzzahl mit "Etage" konvertieren

       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

Es wurde so ein Wert.

Wenn Sie dies jedoch tun möchten, können Sie NumPy verwenden. Es ist schwer zu sagen, dass es für die stochastische Modellierung nicht mehr richtig ist.

Gibt es eine gute Möglichkeit, dieses Problem mit Stan zu lösen?

Zusammenfassung

Teilweise konnten wir Stan verwenden, um das Risiko des neuen Koronavirus abzuschätzen.

Aus der Situation der Diamantprinzessin geht hervor, dass COVID-19 in Bezug auf die Letalität bösartiger ist als die traditionelle Influenza. Eine weitere Bedrohung besteht darin, dass viele Menschen bereits Antikörper gegen Influenza haben, aber nur wenige COVID-19-Antikörper. Sie sollten wachsam sein für die Gesellschaft als Ganzes und nicht für den Einzelnen.

Es wird jedoch geschätzt, dass etwa 1,67 Menschen sterben werden, wenn sie in dieser Größenordnung mit Influenza infiziert sind. Wenn Sie also Angst vor COVID-19 haben, sollten Sie bei normaler Influenza vorsichtiger sein. Aus meiner persönlichen Sicht möchte ich eine Erklärung abgeben.

Fortsetzung stellt auch die Aufgabe in Frage, die Letalität von COVID-19 zu ermitteln. Bitte schau es dir an.

Blinddarm

Quellcode: 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

Ist die neue Corona wirklich eine Bedrohung? Validiert mit Stan (war)
[Fortsetzung] Ist die neue Corona wirklich eine Bedrohung? Berechnete Letalität mit Stan
[New Corona] Ist der nächste Höhepunkt im Dezember? Ich habe die Trendanalyse mit Python versucht!
Erstellen Sie eine neue CSV mit Pandas basierend auf der lokalen CSV
Das Bild ist Namekuji
Ich habe versucht, die Anzahl der im Inland infizierten Menschen der neuen Korona mit einem mathematischen Modell vorherzusagen
Versuchen Sie, mit einem Baseball-Simulator Nr. 1 zu spielen. Ist der Feed Bunt effektiv?
Fehler mit pip: Beim Bestätigen des SSL-Zertifikats ist ein Problem aufgetreten
[Python] Was ist eine with-Anweisung?
Mit PyEphem ist der Weltraum gefährlich
Erstellt eine neue Corona-App in Kyoto mit Pythons Webframework Dash
Das LXC Web Panel, das LXC mit einem Browser bedienen kann, war wunderbar
Ich habe versucht, mit einem Foto einfach ein hochpräzises 3D-Bild zu erstellen [-1]. (Ist der versteckte Bereich wirklich sichtbar?)