python求曲线拐点_使用Python检测新冠肺炎疫情拐点
参考链接：使用Python检测新冠肺炎疫情拐点，抗疫成果明显
1 简介
拐折检测（ Knee fold detection ），即为对于具有上升或下降趋势的曲线而言，在达到该转折点后整体趋势发生明显转变时所存在的那个关键转折位置被称为拐折点（例如，在图1中标记为蓝色的部分之后曲线急剧上升）
图1
本文旨在介绍Python中用于change point detection的第三方包kneed，并分析新型冠状肺炎cases中的各个指标在mathematical sense上的拐点。
2 基于kneed的拐点检测
许多算法普遍采用肘部法则来确定关键参数，在例如在机器学习领域中被广泛应用的聚类分析中，则有明确的应用场景和指导原则。具体而言，在`K-means算法中使用聚类数目k（即簇的数量），而在DBSCAN算法中采用搜索半径eps$作为核心参数进行密度聚类分析。
在面对需明确判断肘部即拐点时 采用基于数学原理的方法进行检测更为科学化
researchers such as Jeannie Albrecht et al. 在 GitHub 仓库开源的论文 Finding a “Kneedle” in a Haystack: Detecting Knee Points in System Behavior 中基于曲率的思想提出了一个系统性地提出了一种拐点检测方案，并将其应用于多种算法研究中。该方案能够有效结合离散与连续数据的特点，在实际应用中展现出良好的效果，并且能够很好地适应不同场景下的数据处理需求。
kneed就是对这篇论文所提出算法的实现。
使用pip install kneed完成安装之后，下面我们来了解其主要用法：
2.1.1 KneeLocator
KneeLocator是kneed中用于检测拐点的模块，其主要参数如下：
x：待检测数据对应的横轴数据序列，如时间点、日期等
y：待检测数据序列，在x条件下对应的值，如x为星期一，对应的y为降水量
S：float型，默认为1，敏感度参数，越小对应拐点被检测出得越快
curve：str型，指明曲线之上区域是凸集还是凹集，concave代表凹，convex代表凸
direction：str型，指明曲线初始趋势是增还是减，increasing表示增，decreasing表示减
online：bool型，用于设置在线/离线识别模式，True表示在线，False表示离线；在线模式下会沿着x轴从右向左识别出每一个局部拐点，并在其中选择最优的拐点；离线模式下会返回从右向左检测到的第一个局部拐点KneeLocator在接收输入参数并完成计算后，在程序中主要关注并返回的属性有哪些？
knee及elbow：返回检测到的最优拐点对应的x
knee_y及elbow_y：返回检测到的最优拐点对应的y
all_elbows及all_knees：返回检测到的所有局部拐点对应的x
all_elbows_y及all_knees_y：返回检测到的所有局部拐点对应的ycurve与direction参数非常重要，用它们组合出想要识别出的拐点模式。
以余弦函数为例，在online设为True时，在以下四种参数组别下对余弦曲线进行拐点计算：曲线=凹形并方向=递增、曲线=凹形并方向=递减、曲线=凸形并方向=递增以及曲线=凸形并方向=递减
import matplotlib.pyplot as plt
from matplotlib import style
import numpy as np
from kneed import KneeLocator
style.use('seaborn-whitegrid')
x = np.arange(1, 3, 0.01)*np.pi
y = np.cos(x)
# 计算各种参数组合下的拐点
kneedle_cov_inc = KneeLocator(x,
y,
curve='convex',
direction='increasing',
online=True)
kneedle_cov_dec = KneeLocator(x,
y,
curve='convex',
direction='decreasing',
online=True)
kneedle_con_inc = KneeLocator(x,
y,
curve='concave',
direction='increasing',
online=True)
kneedle_con_dec = KneeLocator(x,
y,
curve='concave',
direction='decreasing',
online=True)
fig, axe = plt.subplots(2, 2, figsize=[12, 12])
axe[0, 0].plot(x, y, 'k--')
axe[0, 0].annotate(s='Knee Point', xy=(kneedle_cov_inc.knee+0.2, kneedle_cov_inc.knee_y), fontsize=10)
axe[0, 0].scatter(x=kneedle_cov_inc.knee, y=kneedle_cov_inc.knee_y, c='b', s=200, marker='^', alpha=1)
axe[0, 0].set_title('convex+increasing')
axe[0, 0].fill_between(np.arange(1, 1.5, 0.01)*np.pi, np.cos(np.arange(1, 1.5, 0.01)*np.pi), 1, alpha=0.5, color='red')
axe[0, 0].set_ylim(-1, 1)
axe[0, 1].plot(x, y, 'k--')
axe[0, 1].annotate(s='Knee Point', xy=(kneedle_cov_dec.knee+0.2, kneedle_cov_dec.knee_y), fontsize=10)
axe[0, 1].scatter(x=kneedle_cov_dec.knee, y=kneedle_cov_dec.knee_y, c='b', s=200, marker='^', alpha=1)
axe[0, 1].fill_between(np.arange(2.5, 3, 0.01)*np.pi, np.cos(np.arange(2.5, 3, 0.01)*np.pi), 1, alpha=0.5, color='red')
axe[0, 1].set_title('convex+decreasing')
axe[0, 1].set_ylim(-1, 1)
axe[1, 0].plot(x, y, 'k--')
axe[1, 0].annotate(s='Knee Point', xy=(kneedle_con_inc.knee+0.2, kneedle_con_inc.knee_y), fontsize=10)
axe[1, 0].scatter(x=kneedle_con_inc.knee, y=kneedle_con_inc.knee_y, c='b', s=200, marker='^', alpha=1)
axe[1, 0].fill_between(np.arange(1.5, 2, 0.01)*np.pi, np.cos(np.arange(1.5, 2, 0.01)*np.pi), 1, alpha=0.5, color='red')
axe[1, 0].set_title('concave+increasing')
axe[1, 0].set_ylim(-1, 1)
axe[1, 1].plot(x, y, 'k--')
axe[1, 1].annotate(s='Knee Point', xy=(kneedle_con_dec.knee+0.2, kneedle_con_dec.knee_y), fontsize=10)
axe[1, 1].scatter(x=kneedle_con_dec.knee, y=kneedle_con_dec.knee_y, c='b', s=200, marker='^', alpha=1)
axe[1, 1].fill_between(np.arange(2, 2.5, 0.01)*np.pi, np.cos(np.arange(2, 2.5, 0.01)*np.pi), 1, alpha=0.5, color='red')
axe[1, 1].set_title('concave+decreasing')
axe[1, 1].set_ylim(-1, 1)
# 导出图像
plt.savefig('图2.png', dpi=300)图中红色区域分别代表符合参数条件的搜索范围；蓝色三角形则为每种参数组合下通过kneed算法检测出的最佳分界点：
图2
下面我们扩大余弦函数中x的范围，绘制出提取到的所有局部拐点：
x = np.arange(0, 6, 0.01)*np.pi
y = np.cos(x)
# 计算convex+increasing参数组合下的拐点
kneedle = KneeLocator(x,
y,
curve='convex',
direction='increasing',
online=True)
fig, axe = plt.subplots(figsize=[8, 4])
axe.plot(x, y, 'k--')
axe.annotate(s='Knee Point', xy=(kneedle.knee+0.2, kneedle.knee_y), fontsize=10)
axe.set_title('convex+increasing')
axe.fill_between(np.arange(1, 1.5, 0.01)*np.pi, np.cos(np.arange(1, 1.5, 0.01)*np.pi), 1, alpha=0.5, color='red')
axe.fill_between(np.arange(3, 3.5, 0.01)*np.pi, np.cos(np.arange(3, 3.5, 0.01)*np.pi), 1, alpha=0.5, color='red')
axe.fill_between(np.arange(5, 5.5, 0.01)*np.pi, np.cos(np.arange(5, 5.5, 0.01)*np.pi), 1, alpha=0.5, color='red')
axe.scatter(x=list(kneedle.all_knees), y=np.cos(list(kneedle.all_knees)), c='b', s=200, marker='^', alpha=1)
axe.set_ylim(-1, 1)
# 导出图像
plt.savefig('图3.png', dpi=300)实验结果如图3所示，在进行拐点检测的过程中，默认情况下落在数据序列两端的拐点会被自动识别为无效数据点
2.2 探索新冠肺炎疫情数据
下一步我们将尝试将kneed应用于新冠肺炎数据中, 以探究各个指标数学意义上的转折点是否已经出现.
原始数据来源于该GitHub项目经过整理和分析的结果。该GitHub平台基于丁香园提供的最新疫情数据分析。实时提供城市级别的疫情发展数据分析。
在文章开头的部分提到了我的GitHub仓库对应的存储位置，在此路径下可下载用于后续分析的数据集。数据集的更新时间为：2020-02-18 22:55:07 ，接下来我们将对以下内容进行详细分析。
为了便于后续操作和数据分析需求，在开始处理前我们将数据源文件命名为 DXYArea.csv 并导入系统；随后我们对文件中的数据结构进行详细分析；其中，在读取数据时我们将 _updateTime_ 列提前解析为时间格式以便后续的数据处理需求；
import pandas as pd
raw = pd.read_csv('DXYArea.csv', parse_dates=['updateTime'])
raw.info()
查看其第一行信息：
raw.loc[0,:]
可以看出，在各相关省份及其市辖区的最新统计数据显示出详细的人口信息数据。
我们旨在监测全国范围内的现有确诊病例数量、每日新增病例数以及治愈率和死亡率随着时间（以天为单位）的变化趋势图是否已出现数学意义上的拐点。鉴于武汉市的数据更新具有高度复杂性和独特性，并且考虑到其特殊地位，在后续分析中仅针对非武汉市地区的数据进行研究。
首先我们对所有市取每天最晚一次更新的数据作为当天正式的记录值：
# 抽取updateTime列中的年、月、日信息分别保存到新列中
raw['year'], raw['month'], raw['day'] = list(zip(*raw['updateTime'].apply(lambda d: (d.year, d.month, d.day))))
# 得到每天每个市最晚一次更新的疫情数据
temp = raw.sort_values(['provinceName', 'cityName', 'year', 'month', 'day', 'updateTime'],
ascending=False,
ignore_index=True).groupby(['provinceName', 'cityName', 'year', 'month', 'day'])
.agg({'province_confirmedCount': 'first',
'province_curedCount': 'first',
'province_deadCount': 'first',
'city_confirmedCount': 'first',
'city_curedCount': 'first',
'city_deadCount': 'first'})
.reset_index(drop=False)
# 查看前5行
temp.head()
有了已经处理成熟的数据后，我们随后将关注全国（除武汉市外）的关键指标变化趋势的分析。
为了更好地关注重点指标并实现数据可视化的目标
# 计算各指标时序结果
# 全国（除武汉市外）累计确诊人数
nationwide_confirmed_count = temp[temp['cityName'] != '武汉'].groupby(['year', 'month', 'day'])
.agg({'city_confirmedCount': 'sum'})
.reset_index(drop=False)
# 全国（除武汉市外）累计治愈人数
nationwide_cured_count = temp[temp['cityName'] != '武汉'].groupby(['year', 'month', 'day'])
.agg({'city_curedCount': 'sum'})
.reset_index(drop=False)
# 全国（除武汉市外）累计死亡人数
nationwide_dead_count = temp[temp['cityName'] != '武汉'].groupby(['year', 'month', 'day'])
.agg({'city_deadCount': 'sum'})
.reset_index(drop=False)
# 全国（除武汉市外）每日新增确诊人数，即为nationwide_confirmed_count的一阶差分
nationwide_confirmed_inc_count = nationwide_confirmed_count['city_confirmedCount'].diff()[1:]
# 全国（除武汉市外）治愈率
nationwide_cured_ratio = nationwide_cured_count['city_curedCount'] / nationwide_confirmed_count['city_confirmedCount']
# 全国（除武汉市外）死亡率
nationwide_died_ratio = nationwide_dead_count['city_deadCount'] / nationwide_confirmed_count['city_confirmedCount']
# 绘图
#解决中文显示问题
plt.rcParams['font.sans-serif'] = ['KaiTi']
plt.rcParams['axes.unicode_minus'] = False
fig, axes = plt.subplots(3, 2, figsize=[12, 18])
axes[0, 0].plot(nationwide_confirmed_count.index, nationwide_confirmed_count['city_confirmedCount'], 'k--')
axes[0, 0].set_title('累计确诊人数', fontsize=20)
axes[0, 0].set_xticks(nationwide_confirmed_count.index)
axes[0, 0].set_xticklabels([f"{nationwide_confirmed_count.loc[i, 'month']}-{nationwide_confirmed_count.loc[i, 'day']}"
for i in nationwide_confirmed_count.index], rotation=60)
axes[0, 1].plot(nationwide_cured_count.index, nationwide_cured_count['city_curedCount'], 'k--')
axes[0, 1].set_title('累计治愈人数', fontsize=20)
axes[0, 1].set_xticks(nationwide_cured_count.index)
axes[0, 1].set_xticklabels([f"{nationwide_cured_count.loc[i, 'month']}-{nationwide_cured_count.loc[i, 'day']}"
for i in nationwide_cured_count.index], rotation=60)
axes[1, 0].plot(nationwide_dead_count.index, nationwide_dead_count['city_deadCount'], 'k--')
axes[1, 0].set_title('累计死亡人数', fontsize=20)
axes[1, 0].set_xticks(nationwide_dead_count.index)
axes[1, 0].set_xticklabels([f"{nationwide_dead_count.loc[i, 'month']}-{nationwide_dead_count.loc[i, 'day']}"
for i in nationwide_dead_count.index], rotation=60)
axes[1, 1].plot(nationwide_confirmed_inc_count.index, nationwide_confirmed_inc_count, 'k--')
axes[1, 1].set_title('每日新增确诊人数', fontsize=20)
axes[1, 1].set_xticks(nationwide_confirmed_inc_count.index)
axes[1, 1].set_xticklabels([f"{nationwide_confirmed_count.loc[i, 'month']}-{nationwide_confirmed_count.loc[i, 'day']}"
for i in nationwide_confirmed_inc_count.index], rotation=60)
axes[2, 0].plot(nationwide_cured_ratio.index, nationwide_cured_ratio, 'k--')
axes[2, 0].set_title('治愈率', fontsize=20)
axes[2, 0].set_xticks(nationwide_cured_ratio.index)
axes[2, 0].set_xticklabels([f"{nationwide_cured_count.loc[i, 'month']}-{nationwide_cured_count.loc[i, 'day']}"
for i in nationwide_cured_ratio.index], rotation=60)
axes[2, 1].plot(nationwide_died_ratio.index, nationwide_died_ratio, 'k--')
axes[2, 1].set_title('死亡率', fontsize=20)
axes[2, 1].set_xticks(nationwide_died_ratio.index)
axes[2, 1].set_xticklabels([f"{nationwide_dead_count.loc[i, 'month']}-{nationwide_dead_count.loc[i, 'day']}"
for i in nationwide_died_ratio.index], rotation=60)
fig.suptitle('全国范围（除武汉外）', fontsize=30)
# 导出图像
plt.savefig('图7.png', dpi=300)
接着就到了检测拐点的时候了。
为了使代码更加简洁, 我们首先创建自定义函数, 用于从KneeLocator中曲线参数（curve）与方向参数（direction）的所有组合识别合法拐点输出值及其对应趋势变化类型；若无, 则返回None:
def knee_point_search(x, y):
# 转为list以支持负号索引
x, y = x.tolist(), y.tolist()
output_knees = []
for curve in ['convex', 'concave']:
for direction in ['increasing', 'decreasing']:
model = KneeLocator(x=x, y=y, curve=curve, direction=direction, online=False)
if model.knee != x[0] and model.knee != x[-1]:
output_knees.append((model.knee, model.knee_y, curve, direction))
if output_knees.__len__() != 0:
print('发现拐点！')
return output_knees
else:
print('未发现拐点！')下面我们对每个指标进行拐点搜索。
先来看看累计确诊数 ，经过程序的搜索，并未发现有效拐点：
接着检测累计治愈数 ，发现了有效拐点：
在曲线图上标记出拐点：
knee_info = knee_point_search(x=nationwide_cured_count.index,
y=nationwide_cured_count['city_curedCount'])
fig, axe = plt.subplots(figsize=[8, 6])
axe.plot(nationwide_cured_count.index, nationwide_cured_count['city_curedCount'], 'k--')
axe.set_title('累计治愈人数', fontsize=20)
axe.set_xticks(nationwide_cured_count.index)
axe.set_xticklabels([f"{nationwide_cured_count.loc[i, 'month']}-{nationwide_cured_count.loc[i, 'day']}"
for i in nationwide_cured_count.index], rotation=60)
for point in knee_info:
axe.scatter(x=point[0], y=point[1], c='b', s=200, marker='^')
axe.annotate(s=f'{point[2]} {point[3]}', xy=(point[0]+1, point[1]), fontsize=14)
# 导出图像
plt.savefig('图10.png', dpi=300)
基于其convex+increasing特性可知，在2月5日起累计治愈人数呈现出显著的增长态势。
再来看看累计死亡人数 ：
图11
绘制出其拐点：
图12
自2月5日起开始同步运行，并行的累计死亡人数与累计治愈人数呈现出加速上升的趋势。
就日新增确诊数而言，则识别出两个关键转折点。尽管这一指标的变化轨迹呈现明显波动性特征，在第一个关键时间点前的增长速度有所减缓，在第二个关键时间点后则出现显著降幅扩大现象。通过分析其参数信息，则能推断出在第一个转折期增长趋缓，在第二个转折期后下降速度加快：这表明除武汉以外各地的抗疫成效已经取得了一定进展
图13
康复率 和 死亡率 同样出现了关键转折点，在这一过程中,康复率 出现了显著提升的关键转折, 随着医疗工作者群体的共同努力, 更为有效的治疗方法推动了康复率的增长。
图14
虽然最近一次的转折点显示出持续增长的趋势，在对比其与治愈率的变动程度时可以观察到这一指标的增长幅度相对较小。然而，在具体数据上发现，在该期间仅增加了很小的比例。
图15
基于上述分析的结果表明，在这场特殊的抗疫战役中截至目前以来 除去武汉市以外 其他地区的防疫成效已有一定成效 然而 仍需持续投入更多资源以保障防疫成果。