基于上海发布公号最新公布的疫情数据，本文对微博@萧瑟Henry 提到的预测模型进行复现。
本硕读的都是统计学，也是在上海念的书
，忍不住自己动手实操一下，不能把知识都还给我们都喜欢的徐老师啊。
1. 读取上海疫情数据
import pandas as pdimport numpy as npdata = pd.DataFrame({'日期': ['3月12日','13日','14日','15日','16日','17日','18日','19日','20日','21日','22日','23日','24日','25日','26日','27日','28日','29日'], '确诊人数': [1,41,4,5,8,57,8,17,24,31,4,4,29,38,45,50,96,326], '无症状人数': [64,128,130,197,150,203,366,492,734,865,977,979,1580,2231,2631,3450,4381,5656]})data['总人数']=data['确诊人数']+data['无症状人数'] data['LN总人数']=np.log(data['总人数'])data['日期序列']=data.index+1data.set_index('日期', inplace=True)data
2. 观察感染人数随着时间推移的变化趋势
import matplotlib.pyplot as pltplt.rcParams['font.sans-serif']=['SimHei']fig,ax1 = plt.subplots(figsize = (10,6))ax2 = ax1.twinx()l1, = ax1.plot(data.index,data["总人数"],color='green', marker='o', linestyle='dashed',linewidth=2, markersize=10)l2, = ax2.plot(data.index,data["LN总人数"],color='red', marker='o', linestyle='dashed',linewidth=2, markersize=10)ax1.set_ylabel("y1,总人数")ax2.set_ylabel("y2,LN总人数")plt.legend(handles=[l1,l2,],labels=['总人数','LN总人数'],loc='best')ax1.set_title('感染者总人数变化趋势')plt.show()
可以看出样本区间内呈现出指数级增长趋势。
3. 构造预测模型
其中 t 为自变量，表示时间，y 为因变量，表示感染人数，α 和 β 为需要预估的模型参数，ε 为其他干扰项。对上面模型进行拟合，只需对 y 取对数，这里选取对数 e 来尝试，这样就转换成了对 lny 拟合线性回归模型即可。
从下面的预测结果可得到，预测模型如下：
'''构建模型'''import statsmodels.api as smtrain=data[:16] x = sm.add_constant(train['日期序列']) y = train['LN总人数'] model = sm.OLS(y, x) result = model.fit() print(result.summary())
使用 12 日-27 日共 16 个样本进行模型训练，得到检验 R 平方值为97.1%，拟合度非常高，参数检验 P 值均为 0，模型系数估计显著性也很高。
4. 用模型进行预测
data.loc[:,'预测值']=round(np.exp(data.loc[:,'日期序列']*0.2470+4.1919))data.loc[:,'95%置信下限']=round(np.exp(data.loc[:,'日期序列']*0.223+3.964))data.loc[:,'95%置信上限']=round(np.exp(data.loc[:,'日期序列']*0.271+4.420))data
plt.figure(figsize = (10, 6))plt.plot(data.index,data['总人数'],label='实际值',color='green', marker='o', linestyle='dashed',linewidth=2, markersize=10)plt.plot(data.index,data['预测值'],label='预测值',color='red', marker='*', linestyle='dashed',linewidth=2, markersize=12)plt.title("预测值与实际值差异")plt.legend()plt.show()
可以看到拟合效果很好。从未参与到模型训练的 28 号和 29 号两天数据来看，原预测模型有点保守，指数级增长趋势还在，通过模型预测明天感染人数在 7k-8k 之间。
test = pd.DataFrame({'日期': ['3月30日','3月31日']})test['日期序列']=test.index+19test.loc[:,'预测值']=round(np.exp(test.loc[:,'日期序列']*0.2470+4.1919))test.loc[:,'95%置信下限']=round(np.exp(test.loc[:,'日期序列']*0.223+3.964))test.loc[:,'95%置信上限']=round(np.exp(test.loc[:,'日期序列']*0.271+4.420))test
这个数据大家权当看着玩吧，模型没有考虑样本量、参数分布以及其他可能被忽略的重要因素等。
居家隔离第 12 天了，这两天也感受到了社区管控和服务水平的明显提升，开始按照楼栋分批做核酸检测、发放时蔬等物资，尽可能减交叉感染风险，在近期感染数据飙升的背后，很多"补救"措施已在悄悄进行。
希望大家困在疫情的日子里保持好心情，保护好自己，配合防疫工作，早日迎来拐点，迎接解封的日子啊。
题图来源：网站Pexels
THE END
声明：文中观点不代表本站立场。本文传送门：https://eyangzhen.com/239398.html
