【大数据】美国新冠肺炎疫情分析——正确版（QDU）
虽然这一版本是正确的, 但值得注意的是, 其中一些细节层面的思路仍然在错误版本中有所体现, 因此建议首先阅读错误版本
大数据
大数据
问题分析
已知统计数据显示新冠肺炎患者在2020年1月28日至9月8日期间内分别记录了每日新增病例数、累计确诊病例数、治愈病例数以及死亡病例数等关键指标的数据信息。研究者旨在研究上述时间段内各类数据的变化趋势，并建立相应的数学模型以实现对相关指标发展情况的预测分析。
解决思路
问题缺陷
上周完成的"美国新冠分析"实验的核心思路在于：通过六次多项式拟合的方法获得"新增数量"的变化趋势；基于已知信息计算得出"现存确诊数量"后，并结合"新增数量""总确诊数量""治愈数量""死亡数量"之间的关系建立模型；经过30次迭代运算获取未来一个月的信息数据。尽管该研究原本旨在预测未来"总确诊数量"的变化趋势，但实验过程中却直接通过"新增数量"的变化规律推导出了"总确诊数量"的数据信息；本质上该模型是对未来一个月内"治愈率""死亡率"的变化情况进行预测而非直接预测未来日期内的感染总量；这种做法与研究初衷存在重大偏差并且具有明显的局限性：该模型的高度依赖性在于'新增数量'这一单一因素；如果'新增数量'预测出现偏差将直接影响最终结果的质量；同时采用多项式拟合方式对'新增数量'进行建模也存在不合理之处导致最终结果无法令人信服
改进方案
本次实验的改进主要涉及以下几个方面：
预测方向由上周基于误判的"治愈"与"死亡"目标调整为"累计确诊病例总数"的追踪。
病情传播机制主要由现有确诊病例数量决定，在此基础之上。
上次实验相比，本次实验仍采用根据关系进行迭代的方式计算。
关联已知信息
线性回归的标准形式：
y = Xw+b
采用最小二乘的方式：
平方误差： \sum_{i=1}^N(y_i-x_i^Tw)^2
保证平方误差最小，令式子的导数为零可得：
\hat{w}=(X^TX)^{-1}X^Ty
定义如下若干变量：
- 治疗_n：代表第n天累计治疗病例总数，则有treatment_n = 治愈率_n \times 现存确诊_{n-1} + 治愈_{n-1}
- 致命病例数_n：指第n天新增致命病例数量，则有fatal_cases_n = 死亡率_n \times 现存确诊_{n-1} + 死亡_{n-1}
- 现存确诊病例数量_n：指第n天现有确诊病例数量，则由线性回归模型及迭代计算得出
- 康复率_n：指第n天累计康复率指标，则有recovery_rate_n = (治疗_新病例数 - 治愈_{新病例}) / 现存确诊病例_{前一时期}
- 致死率_新计算值n：指第n天计算出的新致死比率，则有$mortality_rate_new[n] = (致命病例数_增量 - 死亡_{前一时期}) / 现存确诊病例_{前一时期}"
- 新增确诊病例总量n：指第n天新增的确诊病例总数，则有 新增确诊病例总量n = 现存确诊病例n + 治愈n + 致死n
我构建了一个线性回归模型来考虑到前五天"现有确诊病例数量"如何影响当天"现有确诊病例数量"的变化趋势，并将其表示为： \mathbf{w} = \begin{bmatrix} w_1 \\ w_2 \\ w_3 \\ w_4 \\ w_5 \\ w_6 \end{bmatrix},其中\mathbf{a}_i \cdot \mathbf{w}分别对应于第n-i天的各种系数乘积之和（i=1, 2, ..., 5），而\mathbf{a}_6 = b为方程中的常数项。\
确定第n天的确诊病例数量所遵循的方程式如下所示：
其中该方程由权重系数与过去五天的确诊病例数量及常数项相乘后相加得出。\
构建矩阵X, 其中矩阵X由六项组成, 分别为 [第n−5天的确诊病例数量, 第n−4天的确诊病例数量, 第n−3天的确诊病例数量, 第n−2天的确诊病例数量, 第n−1天的确诊病例数量, 常数项]
迭代思路：
解决方案
经过对训练/测试数据分割比例的多轮优化处理后，在现有数据基础上按照...:...的比例将这些数据划分为训练样本与验证样本两部分，并通过这一划分策略使得模型在验证阶段的表现更加理想。基于此方法获得的机器学习模型用于进行预测任务。
与上周实验相同，采用指数函数拟合治愈率和死亡率。
治愈数量曲线如下：
死亡数量曲线如下：
针对早期数据的处理问题，在初始病患数量从无逐步增长的过程中存在一定不确定性的情况下，则采用搁置的方法进行处理
准备好所需变量后开始迭代。
得到最终预测未来一个月总确诊数量的曲线如下：
结果分析
由于不再依赖"新增数量"的数据，并通过将每日"总确诊数量"与过去五天的数据建立关系来进行预测
预测未来50天的总确诊人数变化：
预测未来100天的总确诊人数变化：
从图中可以看出该模型具有一定的合理性。
总结展望
优势：
相较于上周使用的模型，在本次实验中我们通过结合实际数据构建了新的数学模型。其中'总确诊病例数'实际上受到'现有确诊病例数'的影响，并且同时综合考虑了过去几天内每天的情况对当前的数据产生的影响。经过多方面的验证与计算得出的结果显示该方法具有较高的准确性
本次代码比上周实验更加简短，更加方便理解。
局限性：
未全面 account for 自然环境及科学技术等要素对数量的影响，在预测死亡率与治愈率时仅提及了"科学技术的进步"会导致死亡率持续下降直至稳定以及治愈率不断攀升直至稳定（考虑到科学技术仍有进一步发展的潜力，在此情况下可忽略趋近于稳定的情况）
附录（代码）
import numpy as np
from numpy import *
from scipy import linalg
from scipy.optimize import curve_fit
from xlrd import open_workbook
import matplotlib.pyplot as plt
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
#标准线性回归函数
def standRegres(xArr, yArr):
xMat = mat(xArr)
yMat = mat(yArr).T
xTx = xMat.T * xMat
#判断行列式为零，则无法求逆
if linalg.det(xTx) == 0:
print('the matrix is singular, cannot do inverse')
return
ws = (xTx).I * (xMat.T*yMat)
return ws
# 根据拟合出的系数计算对应的y值
def linearFunc(ws, xArr):
ws = [item[0] for item in ws.tolist()]
return sum(np.array(ws[:5])*np.array(xArr)) + ws[-1:][0]
# 指数函数，用于拟合变化率
def ExponentialFunction(x, a, b, c):
return a * np.exp(-b * x) + c
# 划分数据集的函数
def splitDataset(x, y, test_size = 0.2):
num = len(y)
train_number = round(num*(1-test_size))
X_train = x[0:train_number]
Y_train = y[0:train_number]
X_test = x[train_number:]
Y_test = y[train_number:]
return X_train, X_test, Y_train, Y_test
# 获取数据
original_data = open_workbook(r'C:\Users\23343\Desktop\test.xlsx') # 打开文件
table = original_data.sheets()[0]
col1 = getcolumn(table, 1)
col2 = getcolumn(table, 2)
col3 = getcolumn(table, 3)
col4 = getcolumn(table, 4)
alive = (array(col2) - array(col3) - array(col4)).tolist()
quantity = len(col1)
index = range(quantity)
# 回归获取系数（对现存确诊进行预测）
xArr = [alive[i:i+5] for i in range(quantity-5)] # 每五个划分为一个列表（注意最后一组没有对应的y，因此range(quantity-5)，而不是range(quantity-4)）
xArr = append(xArr, ones((len(xArr), 1)), axis=1).tolist() # 为每个项添加 1
xArr = [[int(i) for i in item] for item in xArr] # 将浮点转为整数
yArr = [alive[i] for i in range(5, quantity)]
X_train, X_test, Y_train, Y_test = splitDataset(xArr, yArr, 0.23) # 划分数据集
ws = standRegres(X_train, Y_train) # 系数
print(ws)
# """
# 绘图显示效果
xMat = mat(xArr)
yMat = mat(yArr)
yHat = xMat*ws
plt.plot(index[5:len(Y_train)+5], Y_train, 'b*', markersize = '2')
plt.plot(index[len(Y_train)+5:], Y_test, 'r*', markersize = '2')
plt.plot(index[5:], yHat)
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
plt.legend(["训练集", "测试集", "线性回归"])
plt.title("训练与测试", fontsize = '16')
plt.show()
# """
# 设置预测的天数，于此修改，下面预测的天数都会修改
predict_days = 30
# -------- 治愈率预测 --------
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]
popt, pcov = curve_fit(ExponentialFunction, idx_notzero, rate_recovery_discretization, bounds=([0, -1, 0], [np.inf, 0, 50])) # 指数函数拟合
predict_recovery_rate = ExponentialFunction(range(len(index), len(index)+predict_days), popt[0], popt[1], popt[2]) # 预测未来一个月的治愈率
# -------- 死亡率预测 --------
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]
popt, pcov = curve_fit(ExponentialFunction, idx_notzero, rate_death_discretization, bounds=([0, 0, 0], [np.inf, 1, 10])) # 指数函数拟合
predict_death_rate = ExponentialFunction(range(len(index), len(index)+predict_days), popt[0], popt[1], popt[2]) # 预测未来一个月的死亡率
# 变量准备
allpatients = col2.copy()
currentpatients = alive.copy()
predict_allpatients = []
predict_currentpatients = []
predict_recovery_rate
predict_death_rate
predict_recovery = col3[-1:]
predict_death = col4[-1:]
# 开始迭代
for i in range(predict_days):
predict_recovery.append(predict_recovery[i] + predict_recovery_rate[i] * alive[-1:][0])
predict_death.append(predict_death[i] + predict_death_rate[i] * alive[-1:][0])
currentpatients.append(linearFunc(ws, currentpatients[-5:]))
predict_currentpatients.append(currentpatients[-1:][0])
allpatients.append(currentpatients[-1:][0] + predict_recovery[-1:][0] + predict_death[-1:][0])
predict_allpatients.append(allpatients[-1:][0])
plt.plot(range(quantity), allpatients[:-predict_days], 'b*', markersize = '2')
plt.plot(range(quantity, quantity+predict_days), predict_allpatients, 'r*', markersize = '2')
plt.legend(["已知样本点", "预测样本点"])
plt.title("预测未来" + str(predict_days) + "天的总确诊人数", fontsize = '16')
plt.show()