【大数据】美国新冠肺炎疫情分析——错误版（QDU）
这个版本可能存在一些问题，在整体框架上可能还存在一些细节上的不足之处。正确的版本可以直接借鉴这个错误版本中的某些具体思路，在完成之后可能会更加完善和严谨。因此，在开始正式作业之前，请确保自己已经充分理解了这个基础性作业的内容和要求（其实当时老师觉得我们做的还不够深入透彻），然后重新梳理一下自己的思考过程（这样才有机会对自己的分析进行更深层次的反思），这样才能保证最终提交的作品具有较高的质量
大数据
大数据
大数据
大数据
问题分析
已收集自2020年1月28日至2020年9月8日新冠肺炎患者的新增病例数、累计确诊病例数、治愈病例数以及死亡病例数等数据，并分析连续几个月内这些指标的变化规律及演变过程。
解决思路
该思考过程被划分为三个逐步深入的阶段：第一阶段采用多项式拟合方法并进行单独预测；第二阶段则着重分析数量变化率而非具体数值；第三阶段通过整合四种已知信息来构建最终模型。
多项式拟合独立分析
最初未充分考虑四种信息间的关联关系。仅对单个信息进行了多项式拟合以预测趋势。这种思路的优势在于模型结构清晰且独立于其他因素。但存在与实际数据严重偏离的问题，并导致预测结果出现明显偏差。
分析数量变化率而非数量变化
针对“治愈率”和”死亡率“进行预测具有一定的合理性：
- 治愈数量与死亡数量的变动均可通过治愈率与死亡率来体现。
- 变化率与样本总量存在关联性，在对比现有确诊患者时能更好地反映出治愈数量与死亡数量所占的比例关系。
- 对变化率进行分析更加贴近实际情况。
关联四种已知信息并构建最终模型
单独对这四种信息进行分别预测是没有价值的；因为它们之间存在相互关联性的影响因素；所以构建模型时必须基于'新增'数量与'总确诊'数量之间的关系，并结合'总治愈'数量和'总死亡'数量的变化趋势。
在关联四种信息之前，我们需要先了解四者之间的关系。
先定义如下若干变量：
- 治愈_n：代表第n天的确诊治愈总数，
- 死亡_n：代表第n天的确诊死亡总数，
- 现存确诊_n：即为第n天的确诊现存人数，
- 治愈率_n：即第n天的确诊治愈率，
- 死亡率_n：即第n天的确诊死亡率，
- 新增新数_{}: 表示新增加的确诊病例数值，
- 总确诊_{}: 表示累计确诊病例总数，
即满足关系式：
总确诊_{} = 总确诊_{前一时期} + 新增新数
| 日期\类型 | 新增 | 总确诊 | 治愈 | 死亡 |
|---|---|---|---|---|
| ... | ... | ... | ... | ... |
| n-1 | 新增n-1 | 总确诊n-1 | 治愈n-1 | 死亡n-1 |
| n | 新增n | 总确诊n | 治愈n | 死亡n |
假设我们采用某种方法拟合出新增数量曲线, 即未来一个月内新增人数均可预测, 根据上述公式可以通过迭代计算得出未来一个月的总确诊人数。
在第n−1天已知四种信息的情况下，则可通过下述公式计算得出第n−1天的确诊人数：现存确诊_{n−1}=总确诊_{n−1}−治愈_{n−1}−死亡_{n−1}随后我们采用了特定方法来预测未来一个月内的治愈率为R_0与病死率为F_0进而获得每日更新后的相关数据；接着利用上述公式X的变化率为{n}=(X_n − X{n−₁})/现存确诊{n−₁}可迭代计算每日更新后的感染人数变化情况
解决方案
第一阶段：多项式拟合独立预测
基于该思路构建模型，并分别对各类数据应用多项式拟合方法以预测未来一个月的各项指标。
多项式拟合的实现
对给定数据(x_i, y_i)\space\space\space(i=0,1,...,m)，一共m+1个数据点，取多项式P(x)，使
\sum_{i=0}^m r_i^2 = \sum_{i=0}^m [P(x_i)-y_i]^2 = min
函数P(x)称为拟合函数或最小二乘解，设n次多项式P_n(x)=\sum_{k=0}^n a_kx^k，使得
I = \sum_{i=0}^m [P_n(x_i)-y_i]^2 = \sum_{i=0}^m(\sum_{k=0}^n a_kx_i^k-y_i)^2 = min
其中，a_0，a_1，a_2，…，a_n为待求未知数，即多项式每一项的系数，n为多项式的最高次幂，由此，该问题化为求
I = I(a_0, a_1,..., a_n)
的极值问题。由多元函数求极值的必要条件：
\frac{\partial I}{ \partial a_j} = 2\sum_{i=0}^m(\sum_{k=0}^n a_kx_i^k-y_i)x_i^j \space\space\space\space\space\space\space\space j=0, 1, ..., n
变形后得到：
\sum_{k=0}^n(\sum_{i=0}^mx_i^{j+k})a_k = \sum_{i=0}^m x_i^jy_i \space\space\space\space\space\space\space\space j=0, 1, ..., n
这是一个关于a_0，a_1，a_2，…，a_n的线性方程组，用矩阵表示如下：
\begin{bmatrix} m+1 & \sum_{i=0}^mx_i & ... & \sum_{i=0}^m x_i^n \\ \sum_{i=0}^mx_i & \sum_{i=0}^mx_i^2 & ... & \sum_{i=0}^mx_i^{n+1} \\ \vdots & \vdots & & \vdots \\ \sum_{i=0}^mx_i^n & \sum_{i=0}^mx_i^{n+1} & ... & \sum_{i=0}^mx_i^{2n} \\ \end{bmatrix} \space \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_n \\ \end{bmatrix} \space = \space \begin{bmatrix} \sum_{i=0}^my_i \\ \sum_{i=0}^mx_iy_i \\ \vdots \\ \sum_{i=0}^mx_i^my_i \\ \end{bmatrix}
由此可见，在提供数据点(x_i,y_i)以及数量m的基础上，并提供拟合所需的参数n后，则即可确定未知数矩阵(a_0,a_1,...,a_n)。
新增数量的预测曲线
方法一：采用六次多项式进行拟合
在短期内的预测结果较为合理，然而，在更长期的时间范围内进行预测可能会导致计算资源的巨大压力。整体效果欠佳。（对于小于0的点默认为0）
方法二：采用一次多项式拟合（线性回归）
从左侧图表可以看出, 训练表现欠佳; 对于后续的变化趋势表现得较为平稳; 仅能预判大致的变化方向.
总确诊数量的预测曲线
方法一：采用三次多项式拟合
训练集中的样本数量稳定，并且由此使得预测结果相对准确。
方法二：采用逻辑斯蒂回归
在理论状态下病毒传播无限制地进行，并呈指数级增长。然而，在实际情况中增长率难以持续稳定。由于在传播过程中会遇到一定的阻碍，并呈现S型曲线。
某物种迁入新生态系统后其数量会发生变化。若该物种最初数目低于环境的最大承载量，则数目会增加。由于该物种在生态系统中面临天敌、食物、空间等资源不足（非理想条件），则其增长函数符合逻辑斯谛方程……这样的曲线呈S形并被广泛应用于描述资源有限环境中的种群增长模式。而其中的关键公式是：\frac{dx}{dt}=rx(1-x)这个速率参数r决定了生长的速度。而总群体的数量P(t)可以用下面这个公式来表示：P(t) = \frac{KP_0e^{rt}}{K+P_0(e^{rt}-1)}
- K表示环境承载力，在理论生态学中被称为最终能够稳定达到的最大数量。
- $P_0代表初始数量, 即t=0时的状态.
- $r代表种群的相对增长率, 其值越大则种群增长速度越快, 能够更快地向环境承载力靠近; 当其值较小时则种群数量变化较缓慢, 逐渐趋于平衡状态。
该增长模式形成的曲线也被视为s型曲线。在图中左侧展示的是该增长模型下的数量变化情况，在右侧则显示了其增长率趋势。
使用逻辑斯蒂模型对总确诊数量进行预测图如下：
（首先对样本点进行了分段处理，在步长设置为5的情况下实现数据离散化；这种做法避免了过于追求精确度所带来的计算复杂度增加问题，在非线性拟合过程中能够有效减少参数求解难度；随后通过离散化处理后得到较为稀疏的数据点集进而进行拟合运算）
总治愈数量的预测曲线
方法一：采用三次多项式拟合
训练集的样本变化平稳，因此预测结果也比较合理。
总死亡数量的预测曲线
方法一：采用三次多项式拟合
训练集存在波动导致对未来一个月的预测骤减，预测结果较差。
第二阶段：指数函数拟合变化率
总治愈数量的预测曲线
方法二：采用指数函数预测“治愈率”
设指数函数为F(x)=a·e^{-bx}+c，其中a、b和c为参数，通过非线性拟合得到。
总死亡数量的预测曲线
方法二：采用指数函数预测“死亡率”
现实中的治愈比率与死亡比率都应呈现指数变化的趋势，在技术的不断发展过程中，治愈比率将持续提高而死亡比率则将持续下降。
第三阶段：关联四种已知信息并构建最终模型
基于模型分析得出未来每天的新增病例数、治愈率以及死亡率数据。通过持续计算得出未来每天的现存确诊数量、总确诊数量、总治愈数量以及总死亡数量。
最终预测曲线如下：
结果分析
整体模型的预测效果较为理想，在某种程度上表现出较高的合理性。然而，在对'新增数量'这一关键指标进行预判时仍面临较大的挑战。具体而言，在利用六次多项式拟合方法来推算'新增数量'的过程中尽管对未来一个月的结果具有一定的合理性但这种模式未能有效规避可能出现的增长爆炸现象这直接导致后续阶段的预测结果显著下降。此外该指标（即'新增数量'）对于其他相关数据项（如需求量等）的表现具有至关重要的影响其准确度将直接影响整体系统的运行效能
对未来100天的“新增数量”进行六次多项式拟合预测：
从图中可以看出预测100天已经呈爆炸式增长，与现实严重不符。
总结展望
优势
- 经过对变化率的系统性研究与分析
- 利用四种信息之间的相互联系对其进行深度建模
- 采用迭代优化的方法提升模型合理性的效果
局限性
- 未考虑环境等因素对数量的影响
- “新增数量”的预测具有一定的不合理性
展望
- 能够将这些因素纳入模型用于预测
- 能够采用更为先进的模型以实现对'新增数量'的准确预测
附录（代码）
用于调用子函数进行拟合，得到预测的数据后进行迭代并绘制图线。
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from xlrd import open_workbook
from NPolynomialFit import NPolynomialFit, calculate
from getFigure import getFigure
from getcolumn import getcolumn
from logistic import logistic_predict
def ExponentialFunction(x, a, b, c):
return a * np.exp(-b * x) + c
original_data = open_workbook(r'C:\Users\23343\Desktop\test.xlsx') # 打开文件
table = original_data.sheets()[0]
alldates = getcolumn(table, 0) # 获取第一列的日期
col1 = getcolumn(table, 1)
col2 = getcolumn(table, 2)
col3 = getcolumn(table, 3)
col4 = getcolumn(table, 4)
xi = range(len(alldates))
# ------------ 新增 ------------
# ----- 方法一 -----
increase, eachdayinc = getFigure(6, xi, col1, 0, len(alldates)) # k为六次多项式的系数 # 预测后续变化会数量爆炸，即预测效果不佳
# ----- 方法二 -----
# getFigure(1, xi, col1, 0, len(alldates)) # k为一次多项式的系数 # 训练效果不佳，预测后续变化相对平稳，即预测效果还可以
# ------------ 总确诊 ------------
# ----- 方法一 -----
# getFigure(3, xi, col2, 0, len(alldates)) # k为三次多项式的系数 # 预测后续变化相对合理，即预测效果还可以
# ----- 方法二 -----
# xi = [int(xi[i]/5) for i in range(0, len(xi), 5)] # 取点离散化，过于精确反而不利于模型训练与预测
# col2 = [col2[i] for i in range(0, len(col2), 5)]
# logistic_predict(xi, col2) # 逻辑斯蒂回归
# ----- 方法三 -----
# 对“新增”使用六次多项式拟合后，根据新增计算总确诊
# 这种方法的好处在于将“新增”与“总确诊”的关系加入模型中
# plt.figure(2)
# plt.plot(xi, col2, 's', markersize='2') # 原始样本点
# plt.plot(range(len(xi), len(xi)+30), np.array(col2[len(col2)-1]) + increase, 'sr', markersize='2') # 预测未来一个月的总确诊
# plt.legend(["原始样本点", "预测"])
# plt.title("原始数据与预测", fontsize = '16')
# plt.show()
predict_allpatients = np.array(col2[len(col2)-1]) + increase
# ------------ 计算现存患者数 ------------
alive = np.array(col2) - np.array(col3) - np.array(col4)
# ------------ 总治愈 ------------
# ----- 方法一 -----
# getFigure(3, xi, col3, 0, len(alldates)) # k为三次多项式的系数 # 预测后续变化相对合理，即预测效果还可以
# ----- 方法二 -----
# 指数函数拟合，考虑到了实际情况，死亡率同
# 重点看治愈率，对治愈率进行预测
# plt.figure(2)
rate_recovery = [(col3[i + 1] - col3[i]) / alive[i] for i in range(len(col3)-1)]
rate_recovery = np.array(rate_recovery) # 治愈率 = (当日治愈总数 - 前日治愈总数)/前日现存确诊人数
idx_notzero = np.where((rate_recovery != 0) & (rate_recovery < 0.075))[0] # 去掉前面误差较大的点（最开始治愈率低是因为未找到治愈方式）
idx_notzero = [idx_notzero[i] for i in range(0, len(idx_notzero))]
rate_recovery_discretization = [rate_recovery[i] for i in idx_notzero]
# plt.plot(range(1, len(xi)), rate_recovery, 'sr', markersize = '2')
popt, pcov = curve_fit(ExponentialFunction, idx_notzero, rate_recovery_discretization, bounds=([0, -1, 0], [np.inf, 0, 50])) # 指数函数拟合
# plt.plot(range(1, len(xi)+30), ExponentialFunction(range(1, len(xi)+30), popt[0], popt[1], popt[2]), '--')
# plt.plot(range(len(xi), len(xi)+30), ExponentialFunction(range(len(xi), len(xi)+30), popt[0], popt[1], popt[2]), 'sb', markersize = '2')
predict_recovery_rate = ExponentialFunction(range(len(xi), len(xi)+30), popt[0], popt[1], popt[2]) # 预测未来一个月的治愈率
# plt.legend(["样本点","拟合曲线","预测"])
# plt.title("治愈总数的曲线变化预测")
# plt.show()
# ------------ 总死亡 ------------
# ----- 方法一 -----
# getFigure(3, xi, col4, 0, len(alldates)) # k为三次多项式的系数 # 预测后续变化骤降，即效果不佳
# ----- 方法二 -----
# 重点看死亡率，对死亡率进行预测
rate_death = [(col4[i + 1] - col4[i]) / alive[i] for i in range(len(col4) - 1)]
rate_death = np.array(rate_death) # 治愈率 = (当日治愈总数 - 前日治愈总数)/前日现存确诊人数
idx_notzero = np.where((rate_death != 0))[0] # 去掉前面误差较大的点
idx_notzero = [idx_notzero[i] for i in range(0, len(idx_notzero))]
rate_death_discretization = [rate_death[i] for i in idx_notzero]
# plt.plot(range(1, len(xi)), rate_death, 'sr', markersize ='2')
popt, pcov = curve_fit(ExponentialFunction, idx_notzero, rate_death_discretization, bounds=([0, 0, 0], [np.inf, 1, 10])) # 指数函数拟合
# plt.plot(range(1, len(xi)+30), ExponentialFunction(range(1, len(xi)+30), popt[0], popt[1], popt[2]), '--')
# plt.plot(range(len(xi), len(xi)+30), ExponentialFunction(range(len(xi), len(xi)+30), popt[0], popt[1], popt[2]), 'sb', markersize = '2')
predict_death_rate = ExponentialFunction(range(len(xi), len(xi)+30), popt[0], popt[1], popt[2]) # 预测未来一个月的死亡率
# plt.legend(["样本点","拟合曲线","预测"])
# plt.title("死亡总数的曲线变化预测")
# plt.show()
# ------------ 用以下四个信息进行迭代 ------------
predict_allpatients = np.array(predict_allpatients)
predict_recovery_rate = np.array(predict_recovery_rate)
predict_death_rate = np.array(predict_death_rate)
alive = [alive[len(alive)-1]] # 第n-1年现存患者数（仅方便自己理解）
recovery = [col3[len(col3)-1]] # 第n-1年的治愈总数（仅方便自己理解）
death = [col4[len(col4)-1]] # 第n-1年的死亡总数（仅方便自己理解）
predict_days = 30
for days in range(predict_days):
recovery.append(predict_recovery_rate[days] * alive[days] + recovery[days])
death.append(predict_death_rate[days] * alive[days] + death[days])
alive.append(predict_allpatients[days] - recovery[days+1] - death[days+1])
# ------------ 绘制 ------------
plt.subplot(221)
plt.plot(range(0, len(xi)), col1, 'sr', markersize ='2')
plt.plot(range(len(xi), len(xi)+predict_days), eachdayinc, 'sb', markersize='2')
plt.legend(["样本点", "预测点"])
plt.title("新增数量的曲线变化预测", fontsize = "16")
plt.subplot(222)
plt.plot(range(0, len(xi)), col2, 'sr', markersize ='2')
plt.plot(range(len(xi), len(xi)+predict_days), predict_allpatients, 'sb', markersize='2')
plt.legend(["样本点", "预测点"])
plt.title("确诊总数的曲线变化预测", fontsize = "16")
plt.subplot(223)
plt.plot(range(0, len(xi)), col3, 'sr', markersize ='2')
plt.plot(range(len(xi), len(xi)+predict_days), recovery[1:], 'sb', markersize='2')
plt.legend(["样本点", "预测点"])
plt.title("治愈总数的曲线变化预测", fontsize = "16")
plt.subplot(224)
plt.plot(range(0, len(xi)), col4, 'sr', markersize ='2')
plt.plot(range(len(xi), len(xi)+predict_days), death[1:], 'sb', markersize='2')
plt.legend(["样本点", "预测点"])
plt.title("死亡总数的曲线变化预测", fontsize = "16")
plt.show()调用多项式拟合函数并获取预测值进行绘图。
import numpy as np
import matplotlib.pyplot as plt
from NPolynomialFit import NPolynomialFit, calculate
def getFigure(n, xi, yi, l, r, days = 30):
k, x, y, xi_y = NPolynomialFit(n, xi, yi, l, r) # n为多项式拟合次数，xi、yi为样本点，l、r表示预测区间
error = abs(np.array(xi_y) - np.array(yi)) # 计算误差
sub1 = plt.subplot(121)
plt.plot(range(len(xi)+days), np.zeros(len(xi)+days), '--k') # 基准线
plt.plot(xi, yi, 's', markersize='2') # 原始样本点
predict_days = range(len(xi), len(xi)+days)
predict_nums = calculate(k, predict_days)
plt.plot(predict_days, predict_nums, 's', markersize='2')
plt.plot(x, y, '-r') # 拟合曲线
sub1.set_title("原始数据与预测", fontsize = '16')
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
plt.legend(["0-基准线", "原始样本点", "预测", "拟合曲线"])
sub2 = plt.subplot(122)
plt.plot(range(len(xi)+days), np.zeros(len(xi)+days), '--k') # 基准线
plt.plot(xi, error, 'sr', markersize='2')
sub2.set_title("误差分析", fontsize = '16')
plt.show()
return np.array(predict_nums).cumsum(), predict_nums实现n次多项式拟合。
import math
import numpy as np
from numpy.linalg import inv
def calculate(k, x):
y = []
for xi in x:
mi = 0
tmp = 0
for ki in k:
tmp = tmp + ki*math.pow(xi, mi)
mi = mi + 1
y.append(tmp)
return y
def NPolynomialFit(n, xi, yi, l, r):
"""xi，yi表示样本数据集，l，r表示要拟合的区间"""
m = len(xi)-1 # 总共m+1个样本点
A = [[0 for j in range(n+1)] for i in range(n + 1)]
B = [[0] for i in range(n+1)]
for row in range(n + 1):
for col in range(n + 1):
for i in range(m + 1):
A[row][col] = A[row][col] + math.pow(xi[i], row + col) # 多项式系数的“系数”
for row in range(n + 1):
for i in range(m + 1):
B[row][0] = B[row][0] + yi[i] * math.pow(xi[i], row) # 得到等式右侧式子
A = np.array(A)
B = np.array(B)
k = np.dot(inv(A), B) # 计算多项式的系数
k = [i[0] for i in k.tolist()] # 化二维为一维
x = [l + i * 0.001 for i in range((r-l)*1000)] # 以0.001为步长从l取到r
y = calculate(k, x) # 计算对应的y值
xi_y = calculate(k, xi)
return k, x, y, xi_y逻辑斯蒂模型预测并绘图。
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
def logistic_increase_function(t,K,P0,r):
# t:time t0:initial time P0:initial_value K:capacity r:increase_rate
t0 = 0
exp_value=np.exp(r*(t-t0))
return (K*exp_value*P0)/(K+(exp_value-1)*P0)
def logistic_predict(xi, yi):
step=5
xi = np.array(xi)
yi = np.array(yi)
# 调参
popt, pcov = curve_fit(logistic_increase_function, xi, yi, bounds=([35000000, -np.inf, -np.inf], [np.inf, np.inf, np.inf]))
xi_y = logistic_increase_function(xi, popt[0], popt[1], popt[2])
tomorrow = [xi[xi.shape[0]-1] + i for i in range(1, 60)]
tomorrow = np.array(tomorrow)
tomorrow_predict = logistic_increase_function(tomorrow, popt[0], popt[1], popt[2])
# 绘图
plot1 = plt.plot(xi*step, yi, 's', label="confimed infected people number")
plot2 = plt.plot(xi*step, xi_y, 'r', label='predict infected people number')
plot3 = plt.plot(tomorrow*step, tomorrow_predict, 'b', label='predict infected people number')
plt.xlabel('time')
plt.ylabel('confimed infected people number')
plt.legend(loc=0) # 指定legend的位置右下角
print(logistic_increase_function(np.array(28), popt[0], popt[1], popt[2]))
print(logistic_increase_function(np.array(29), popt[0], popt[1], popt[2]))
plt.show()获取每列的数据信息。
from datetime import datetime
def getcolumn(table, i):
res = []
for val in table.col_values(i):
if i == 0:
if len(val) != 2:
date = val.split('.')
res.append(datetime(int(date[0]), int(date[1]), int(date[2])))
else:
try:
res.append(int(val))
except:
continue
return res