项目介绍
这个项目旨在通过利用传染病模型，结合实际观测数据，实现对传染病传播过程的更准确预测。我们采用了多种经典传染病模型，包括SIR、SIR模型带有随机性、SEIR、SEIR模型带有随机性、SI、SIRS、SEIRS-V以及SIRD，并通过优化算法对模型参数进行调整，以最好地拟合现实世界的数据。
关键词
SIR、SIR model with stochasticity、SEIR、SEIR model with stochasticity、SI、SIRS、SEIRS-V、SIRD
项目思路
整个项目的实现逻辑可以分为以下几个关键步骤：
- 数据加载与处理：
- 从CSV文件中加载实际观测数据，包括累积确诊、治愈和死亡数据。
- 进行数据预处理，如索引设置和数据类型转换。
- 传染病模型定义：
- 实现多种传染病模型，包括SIR、SIR模型带有随机性、SEIR、SEIR模型带有随机性、SI、SIRS、SEIRS-V以及SIRD等。
- 每个模型都有相应的微分方程描述传播过程。
- 时间序列扩展：
- 通过扩展时间序列，使其覆盖预测范围，以匹配模型预测的时间点。
- 优化参数：
- 定义损失函数，用于衡量模型预测与实际数据的差异。
- 利用数学优化算法（例如L-BFGS-B）调整模型参数，以最小化损失函数。
- 模型预测：
- 使用优化后的参数，结合传染病模型和ODE求解器，对未来一定时间范围内的传播过程进行预测。
- 得到预测结果，包括感染者、康复者和死亡者的数量。
- 结果可视化：
- 将实际数据和模型预测结果可视化，以直观展示模型拟合效果。
- 使用Matplotlib等工具创建图表，包括感染者、康复者和死亡者随时间的变化趋势。
具体代码
SIR
#!/usr/bin/python
import numpy as np
import pandas as pd
from scipy.integrate import solve_ivp
from scipy.integrate import odeint
from scipy.optimize import minimize
import matplotlib.pyplot as plt
from datetime import timedelta, datetime
predict_range = 200
s_0=99998
i_0=2
r_0=0
class Learner(object):
def __init__(self, loss, predict_range, s_0, i_0, r_0):
self.loss = loss
self.predict_range = predict_range
self.s_0 = s_0
self.i_0 = i_0
self.r_0 = r_0
def load_confirmed(self):
df = pd.read_csv('02_SZ_DailyCases.csv')
df.set_index(["Date"], inplace=True)
dff=df["cummulative confirmed cases"]
print(dff.T)
return dff.T
def load_recovered(self):
df = pd.read_csv('02_SZ_DailyCases.csv')
df.set_index(["Date"], inplace=True)
dff = df["cummulative cured cases"]
print(dff.T)
return dff.T
def load_dead(self):
df = pd.read_csv('02_SZ_DailyCases.csv')
df.set_index(["Date"], inplace=True)
dff = df["cummulative dead cases"]
print(dff.T)
return dff.T
def extend_index(self, index, new_size):
values = index.values
current = datetime.strptime(index[-1], '%m/%d/%y')
while len(values) < new_size:
current = current + timedelta(days=1)
values = np.append(values, datetime.strftime(current, '%m/%d/%y'))
return values
def predict(self, beta, gamma, data, recovered, death, s_0, i_0, r_0):
new_index = self.extend_index(data.index, self.predict_range)
size = len(new_index)
def SIR(t, y):
S = y[0]
I = y[1]
R = y[2]
return [-beta * S * I, beta * S * I - gamma * I, gamma * I]
extended_actual = np.concatenate((data.values, [None] * (size - len(data.values))))
extended_recovered = np.concatenate((recovered.values, [None] * (size - len(recovered.values))))
extended_death = np.concatenate((death.values, [None] * (size - len(death.values))))
return new_index, extended_actual, extended_recovered, extended_death, solve_ivp(SIR, [0, size],
[s_0, i_0, r_0],
t_eval=np.arange(0, size, 1))
def train(self):
recovered = self.load_recovered()
death = self.load_dead()
data = (self.load_confirmed() - recovered - death)
optimal = minimize(loss, [0.001, 0.001], args=(data, recovered, self.s_0, self.i_0, self.r_0),
method='L-BFGS-B', bounds=[(0.00000001, 0.4), (0.00000001, 0.4)])
print(optimal)
beta, gamma = optimal.x
new_index, extended_actual, extended_recovered, extended_death, prediction = self.predict(beta, gamma, data,
recovered, death,
self.s_0, self.i_0,
self.r_0)
df = pd.DataFrame(
{
'Infected data': extended_actual, 'Recovered data': extended_recovered, 'Death data': extended_death,
'Susceptible': prediction.y[0], 'Infected': prediction.y[1], 'Recovered': prediction.y[2]},
index=new_index)
fig, ax = plt.subplots(figsize=(15, 10))
# df.to_csv("result_SIR.csv")
#ax.set_title(self.country)
df.plot(ax=ax)
print(f" beta={
beta:.8f}, gamma={
gamma:.8f}, r_0:{
(beta / gamma):.8f}")
# fig.savefig("result_SIR.png")
def loss(point, data, recovered, s_0, i_0, r_0):
size = len(data)
beta, gamma = point
def SIR(t, y):
S = y[0]
I = y[1]
R = y[2]
return [-beta * S * I, beta * S * I - gamma * I, gamma * I]
# Here
# solution = odeint(SIR, [0, size], [s_0, i_0, r_0])
solution = solve_ivp(SIR, [0, size], [s_0, i_0, r_0], t_eval=np.arange(0, size, 1), vectorized=True)
l1 = np.sqrt(np.mean((solution.y[1] - data) ** 2))
l2 = np.sqrt(np.mean((solution.y[2] - recovered) ** 2))
alpha = 0.1
return alpha * l1 + (1 - alpha) * l2
def main():
learner = Learner(loss, predict_range, s_0, i_0, r_0)
learner.train()
if __name__ == '__main__':
main()
SIR model with stochasticity
#!/usr/bin/python
import numpy as np
import pandas as pd
from scipy.integrate import solve_ivp
from scipy.integrate import odeint
from scipy.optimize import minimize
import matplotlib.pyplot as plt
from datetime import timedelta, datetime
np.random.seed(42) # for reproducibility
predict_range = 200
s_0=99998
i_0=2
r_0=0
class Learner(object):
def __init__(self, loss, predict_range, s_0, i_0, r_0):
self.loss = loss
self.predict_range = predict_range
self.s_0 = s_0
self.i_0 = i_0
self.r_0 = r_0
def load_confirmed(self):
df = pd.read_csv('02_SZ_DailyCases.csv')
df.set_index(["Date"], inplace=True)
dff=df["cummulative confirmed cases"]
print(dff.T)
return dff.T
def load_recovered(self):
df = pd.read_csv('02_SZ_DailyCases.csv')
df.set_index(["Date"], inplace=True)
dff = df["cummulative cured cases"]
print(dff.T)
return dff.T
def load_dead(self):
df = pd.read_csv('02_SZ_DailyCases.csv')
df.set_index(["Date"], inplace=True)
dff = df["cummulative dead cases"]
print(dff.T)
return dff.T
def extend_index(self, index, new_size):
values = index.values
current = datetime.strptime(index[-1], '%m/%d/%y')
while len(values) < new_size:
current = current + timedelta(days=1)
values = np.append(values, datetime.strftime(current, '%m/%d/%y'))
return values
def predict(self, beta, gamma, data, recovered, death, s_0, i_0, r_0, sigma=0.1):
new_index = self.extend_index(data.index, self.predict_range)
size = len(new_index)
def SIR_stochastic
文章评论