-
Notifications
You must be signed in to change notification settings - Fork 0
/
run_infection_daily.py
138 lines (114 loc) · 6.12 KB
/
run_infection_daily.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
import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt
from model import load_from_PC_increment_daily, load_from_PC_delta_infected_daily
import matplotlib.ticker as ticker
import datetime as dt
from datetime import date
from sklearn.metrics import mean_squared_error
DATE_TO_START = '2020-02-24'
DAYS = 80
def main():
t = np.arange(0, DAYS)
# ----------------------------#
# ----- INFECTION DAILY ----- #
infection = load_from_PC_increment_daily(remote=True)
infection_extended = np.concatenate((infection.values, [None] * (DAYS - len(infection.values))))
infection_roll = infection.rolling(5).mean()
infection_roll_extended = np.concatenate((infection_roll.values, [None] * (DAYS - len(infection_roll.values))))
x = np.arange(0, len(infection))
# LOGISTIC FIT
mean_est = sum(x * infection.values) / sum(infection.values)
(a, c, d), _ = opt.curve_fit(logistic_function_pdf, x, infection.values, p0=[173084, 0.13, mean_est])
print(a, c, d)
infection_fit_logistic = logistic_function_pdf(t, a, c, d)
print('MSE: {:.0f}'.format(mean_squared_error(infection_fit_logistic[:len(infection.values)], infection.values)))
# GAUSSIAN FIT
'''
mean_est = sum(x * infection.values) / sum(infection.values)
sigma_est = sum(infection.values * (x - mean_est) ** 2) / sum(infection.values)
(amplitude, mu, std), _ = opt.curve_fit(gaussian_function_pdf, x, infection.values, p0=[max(infection.values), mean_est, sigma_est])
print(amplitude, mu, std)
infection_fit_normal = gaussian_function_pdf(t, amplitude, mu, std)
print('MSE: {:.0f}'.format(mean_squared_error(infection_fit_normal[:len(infection.values)], infection.values)))
'''
# GOMPERTZ FIT
(a,b,n), pcov = opt.curve_fit(shifted_gompertz_function_pdf, x, infection.values, p0=[1000, 1, 0.1])
print(a,b,n)
infection_fit_gompertz = shifted_gompertz_function_pdf(t, a, b, n)
print('MSE: {:.0f}'.format(mean_squared_error(infection_fit_gompertz[:len(infection.values)], infection.values)))
# --------------------------- #
# -- DELTA INFECTION DAILY -- #
delta_infection = load_from_PC_delta_infected_daily(remote=True)
infection_delta_extended = np.concatenate((delta_infection.values, [None] * (DAYS - len(delta_infection.values))))
infection_delta_roll = delta_infection.rolling(5).mean()
infection_delta_roll_extended = np.concatenate((infection_delta_roll.values, [None] * (DAYS - len(infection_delta_roll.values))))
x = np.arange(0, len(delta_infection))
# LOGISTIC FIT
mean_est = sum(x * delta_infection.values) / sum(delta_infection.values)
(a, c, d), pcov = opt.curve_fit(logistic_function_pdf, x, delta_infection.values, p0=[173084, 0.13, mean_est])
print(a, c, d)
infection_delta_fit_logistic = logistic_function_pdf(t, a, c, d)
print('MSE: {:.0f}'.format(mean_squared_error(infection_delta_fit_logistic[:len(delta_infection.values)], delta_infection.values)))
# GAUSSIAN FIT
'''
mean_est = sum(x * delta_infection.values) / sum(delta_infection.values)
sigma_est = sum(delta_infection.values * (x - mean_est) ** 2) / sum(delta_infection.values)
(amplitude, mu, std), pcov = opt.curve_fit(gaussian_function_pdf, x, delta_infection.values, p0=[max(delta_infection.values), mean_est, sigma_est])
print(amplitude, mu, std)
infection_delta_fit_normal = gaussian_function_pdf(t, amplitude, mu, std)
print('MSE: {:.0f}'.format(mean_squared_error(infection_delta_fit_normal[:len(delta_infection.values)], delta_infection.values)))
'''
# GOMPERTZ FIT
(a,b,n), pcov = opt.curve_fit(shifted_gompertz_function_pdf, x, delta_infection.values, p0=[114691, 0.09, 10])
print(a,b,n)
infection_delta_fit_gompertz = shifted_gompertz_function_pdf(t, a, b, n)
print('MSE: {:.0f}'.format(mean_squared_error(infection_delta_fit_gompertz[:len(delta_infection.values)], delta_infection.values)))
# --------------------------- #
# ------- PLOT RESULTS ------ #
def todate(x, pos):
return (dt.datetime.strptime(DATE_TO_START, '%Y-%m-%d') + dt.timedelta(days=x)).strftime('%d/%m')
fmt = ticker.FuncFormatter(todate)
# subplot 1
ax = plt.subplot(211)
ax.plot(t, infection_extended, 'ro', alpha=1, mfc='none', label='Infected')
ax.plot(t, infection_roll_extended,'r--', alpha=0.5, label='Infected MA 5d')
ax.plot(t, infection_fit_logistic, label='Logistic fit')
ax.plot(t, infection_fit_gompertz, label='Gompertz fit')
ax.set_ylabel('Infected Daily')
legend = ax.legend()
legend.get_frame().set_alpha(0.5)
ax.xaxis.set_major_formatter(fmt)
ax.xaxis.set_tick_params(rotation=30)
ax.xaxis.set_major_locator(ticker.MultipleLocator(5))
ax.yaxis.set_tick_params(length=0)
ax.xaxis.set_tick_params(length=0)
ax.grid(b=True)
# subplot 2
ax2 = plt.subplot(212)
ax2.plot(t, infection_delta_extended, 'ro', alpha=1, mfc='none', label='Infected')
ax2.plot(t, infection_delta_roll_extended,'r--', alpha=0.5, label='Infected MA 5d')
ax2.plot(t, infection_delta_fit_logistic, label='Logistic fit')
#ax2.plot(t, infection_delta_fit_normal, label='Normal fit')
ax2.plot(t, infection_delta_fit_gompertz, label='Gompertz fit')
ax2.set_ylabel('Infected Delta Daily')
legend = ax2.legend()
legend.get_frame().set_alpha(0.5)
ax2.xaxis.set_major_formatter(fmt)
ax2.xaxis.set_tick_params(rotation=45)
ax2.xaxis.set_major_locator(ticker.MultipleLocator(5))
ax2.yaxis.set_tick_params(length=0)
ax2.xaxis.set_tick_params(length=0)
ax2.grid(b=True)
plt.savefig('forecast_daily/COVID_forecast_{}'.format(date.today().strftime('%d_%m')))
plt.show()
def logistic_function_cdf(x, a, b, c, d):
return a / (1. + np.exp(-c * (x - d))) + b
def logistic_function_pdf(x, a, c, d):
return a * c * np.exp(-c * (x - d)) / (1. + np.exp(-c * (x - d)))**2
def gaussian_function_pdf(x, amplitude, mean, stddev):
return amplitude * np.exp(-((x - mean) / 4 / stddev)**2)
def shifted_gompertz_function_pdf(x, amplitude, b, n):
return amplitude * b * np.exp(-b * x - n * np.exp(-b*x)) * (1 + n * (1 - np.exp(-b * x)))
if __name__ == '__main__':
main()