he_et_al_regraph.py

import scipy.stats
import math

# From He et al https://www.nature.com/articles/s41591-020-0869-5#MOESM1
# Infection-to-symptom-onset (incubation period) is log-normal distributed with
incubation_period = scipy.stats.lognorm(0.6612, scale=math.exp(1.434065))

import matplotlib.pyplot as plt
import numpy as np
x = np.linspace(incubation_period.ppf(0.001), incubation_period.ppf(0.999), 100)
plt.plot(x, incubation_period.pdf(x), 'k-', lw=2, label='frozen pdf')

# Infectiousness per day after symptom onset is gamma distributed with
infections_prop_by_day_post_symptoms = scipy.stats.gamma(2.1157790, scale=1/0.6898583, loc=-2.3066912)
x = np.linspace(infections_prop_by_day_post_symptoms.ppf(0.001), infections_prop_by_day_post_symptoms.ppf(0.999), 100)
plt.plot(x, infections_prop_by_day_post_symptoms.pdf(x), 'k-', lw=2, label='frozen pdf')

# Serial interval (days from infector's symptom onset to infectee's symptom onset) is gamma distributed with
serial_interval = scipy.stats.gamma(8.1237578, scale=1/0.6361684, loc=-6.5)
x = np.linspace(serial_interval.ppf(0.001), serial_interval.ppf(0.999), 100)
plt.plot(x, serial_interval.pdf(x), 'k-', lw=2, label='frozen pdf')

x = np.linspace(-5, 20, 100)

infectiousness_sample = infections_prop_by_day_post_symptoms.rvs(1000000)
incubation_sample = incubation_period.rvs(1000000)
empirical_serial_interval = infectiousness_sample + incubation_sample
hist, edges = np.histogram(empirical_serial_interval, 25, range=(-5, 20), density=True)

plt.figure(figsize=(7, 7))
plt.subplot(211, xlabel="Days after infection", ylabel="Probability density")
plt.minorticks_on()
plt.plot(x, incubation_period.pdf(x), 'r-', label="Incubation period")
plt.grid(True, axis="x", which="both")
plt.grid(True, axis="y", which="major")
plt.legend()
plt.subplot(212, xlabel="Days after symptom onset", ylabel="Probability density")
plt.minorticks_on()
plt.plot(x, infections_prop_by_day_post_symptoms.pdf(x), 'k-', label="Infectiousness")
plt.plot(x, serial_interval.pdf(x), 'b-', label="Serial interval")
plt.plot(edges[:-1], hist, "b--", label="Implied serial interval")
plt.grid(True, axis="x", which="both")
plt.grid(True, axis="y", which="major")
plt.legend()