-
Notifications
You must be signed in to change notification settings - Fork 33
/
distributions.py
204 lines (164 loc) · 7.94 KB
/
distributions.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
import time
import bisect
import numpy as np
import pandas as pd
import networkx as nx
import scipy
import scipy.optimize
import scipy as sp
import random as rd
import os
import math
import matplotlib
import matplotlib.pyplot as plt
TO_HOURS = 24.0
class CovidDistributions(object):
"""
Class to sample from specific distributions for SARS COV2
"""
def __init__(self, country):
'''
Covid-19 specific constants from literature
ALL UNITS IN DAYS
'''
self.lambda_0 = 0.0
# https://www.medrxiv.org/content/10.1101/2020.03.03.20029983v1
# https://www.who.int/docs/default-source/coronaviruse/who-china-joint-mission-on-covid-19-final-report.pdf
# https://science.sciencemag.org/content/368/6491/eabb6936/tab-pdf
self.R0 = 2.0 # for seeding
# Proportion of infections that are asymptomatic
# https://www.medrxiv.org/content/10.1101/2020.02.03.20020248v2
# https://science.sciencemag.org/content/368/6491/eabb6936/tab-pdf
# https://www.medrxiv.org/content/10.1101/2020.04.17.20053157v1.full.pdf
self.alpha = 0.4
# Relative transmission rate of asymptomatic individuals
# Li et al (Science, 2020): "multiplicative factor reducing the transmission rate of unreported infected patients"
# https://science.sciencemag.org/content/368/6490/489/tab-pdf
self.mu = 0.55
# https://npgeo-corona-npgeo-de.hub.arcgis.com/datasets/dd4580c810204019a7b8eb3e0b329dd6_0
# https://www.imperial.ac.uk/media/imperial-college/medicine/sph/ide/gida-fellowships/Imperial-College-COVID19-NPI-modelling-16-03-2020.pdf
# https://www.bag.admin.ch/bag/en/home/krankheiten/ausbrueche-epidemien-pandemien/aktuelle-ausbrueche-epidemien/novel-cov/situation-schweiz-und-international.html
if country == 'GER':
self.fatality_rates_by_age = np.array([0.0, 0.0, 0.0, 0.004, 0.073, 0.247])
self.p_hospital_by_age = np.array([0.001, 0.002, 0.012, 0.065, 0.205, 0.273])
elif country == 'CH':
self.fatality_rates_by_age = np.array([0, 0, 0, 0.001, 0.001, 0.005, 0.031, 0.111, 0.265])
self.p_hospital_by_age = np.array([0.155, 0.038, 0.028, 0.033, 0.054, 0.089, 0.178, 0.326, 0.29])
else:
raise NotImplementedError('Invalid country requested.')
# https://www.medrxiv.org/content/10.1101/2020.03.09.20033217v2
# self.lab_aerosol_halflife = 2.0
# self.real_site_halflife_factor = 1.0 # 10-times shorter halflife at real sites
self.lab_aerosol_halflife = 1.1 # hours
self.real_site_halflife_factor = 0.1 # 10-times shorter halflife at real sites
self.real_site_aerosol_halflife = self.lab_aerosol_halflife * self.real_site_halflife_factor
# 6.301338005090411
self.gamma = np.log(2.0) / self.real_site_aerosol_halflife
# 0.3654120904376099
self.delta = np.log(1 / 0.10) / self.gamma # time of intensity decrease to below 10 %
# self.delta = np.log(1 / 0.20) / self.gamma
# Incubation period: estimated mean is 5.52 days, std dev is 2.41 days
# To be able to approx. represent latent and infectious incubation period separately,
# we subtract the estimated median infectious time period from the above mean
# when sampling the arrival times. The infectious period std dev is heuristically set to 1 day.
# https://www.medrxiv.org/content/10.1101/2020.03.05.20030502v1
# https://www.nature.com/articles/s41591-020-0869-5
# https://www.who.int/docs/default-source/coronaviruse/who-china-joint-mission-on-covid-19-final-report.pdf
# https://jamanetwork.com/journals/jama/fullarticle/2761044
# https://pubmed.ncbi.nlm.nih.gov/32079150/
self.incubation_mean_of_lognormal = 5.52
self.incubation_std_of_lognormal = 2.41
self.median_infectious_without_symptom = 2.3
# Other literature parameters
self.median_symp_to_resi = 14.0
self.median_asymp_to_resi = 14.0
self.median_symp_to_hosp = 7.0
self.symp_to_death_mean_of_lognormal = 13.0
def normal_to_lognormal(self, mu, std):
'''
Converts mean and std to lognormal mean and std
'''
phi = np.sqrt(np.square(std) + np.square(mu))
lnmu = np.log(np.square(mu) / phi)
lnstd = np.sqrt(np.log(np.square(phi) / np.square(mu)))
return lnmu, lnstd
def __mean_distribution(self, mu, std, size):
'''
Samples from log normal distribution with a specified mean and std dev
'''
lnmean, lnstd = self.normal_to_lognormal(mu, std)
return np.random.lognormal(mean=lnmean, sigma=lnstd, size=size)
def sample_susc_baseexpo(self, size=1):
'''
Samples r.v. of susc -> expo (at base rate only)
'''
# this is base rate exposure only
assert(self.lambda_0 != 0.0)
return TO_HOURS * np.random.exponential(
scale=1.0 / self.lambda_0,
size=size)
def sample_expo_ipre(self, size=1):
'''
Samples r.v. of expo -> ipre
'''
return TO_HOURS * (self.__mean_distribution(
self.incubation_mean_of_lognormal - self.median_infectious_without_symptom,
self.incubation_std_of_lognormal, size=size))
def sample_expo_iasy(self, size=1):
'''
Samples r.v. of expo -> iasy
'''
return TO_HOURS * (self.__mean_distribution(
self.incubation_mean_of_lognormal - self.median_infectious_without_symptom,
self.incubation_std_of_lognormal, size=size))
def sample_ipre_isym(self, size=1):
'''
Samples r.v. of ipre -> isym (Incubation period)
'''
return TO_HOURS * self.__mean_distribution(self.median_infectious_without_symptom, 1.0, size=size)
def sample_isym_resi(self, size=1):
'''
Samples r.v. of isym -> resi
'''
return TO_HOURS * self.__mean_distribution(self.median_symp_to_resi, 1.0, size=size)
def sample_isym_dead(self, size=1):
'''
Samples r.v. of isym -> dead
'''
return TO_HOURS * self.__mean_distribution(self.symp_to_death_mean_of_lognormal, 1.0, size=size)
def sample_isym_hosp(self, size=1):
'''
Samples r.v. of isym -> hosp
'''
return TO_HOURS * self.__mean_distribution(self.median_symp_to_hosp, 1.0, size=size)
def sample_iasy_resi(self, size=1):
'''
Samples r.v. of iasy -> resi
'''
return TO_HOURS * self.__mean_distribution(self.median_asymp_to_resi, 1.0, size=size)
def sample_is_fatal(self, ages, size=1):
'''
Samples iid Bernoulli r.v. of fatality based on age group
'''
assert(ages.shape[0] == size[0])
return np.random.binomial(1, self.fatality_rates_by_age[ages], size=size)
def sample_is_hospitalized(self, ages, size=1):
'''
Samples iid Bernoulli r.v. of hospitalization based on age group
'''
assert(ages.shape[0] == size[0])
return np.random.binomial(1, self.p_hospital_by_age[ages], size=size)
if __name__ == '__main__':
dist = CovidDistributions(country='GER')
print('expo to ipre/iasy (subtracted infectious window before symptoms) : ',
dist.normal_to_lognormal(dist.incubation_mean_of_lognormal - dist.median_infectious_without_symptom, dist.incubation_std_of_lognormal))
print('ipre to isym : ',
dist.normal_to_lognormal(dist.median_infectious_without_symptom, 1.0))
print('isym to resi : ',
dist.normal_to_lognormal(dist.median_symp_to_resi, 1.0))
print('isym to dead : ',
dist.normal_to_lognormal(dist.symp_to_death_mean_of_lognormal, 1.0))
print('isym to hosp : ',
dist.normal_to_lognormal(dist.median_symp_to_hosp, 1.0))
print('iasy to resi : ',
dist.normal_to_lognormal(dist.median_asymp_to_resi, 1.0))