在前一段时间,我脑子一热使用指数回归以及多项式回归对新型冠状病毒2019-nCov的感染人数进行预测。后来发现不行啊,感染人数不可能一直上涨啊,总得有停止上涨的时候啊!在多名网友的提醒下,在参考了邢翔瑞大佬的博客之后,我痛定思痛,尝试使用Logistic增长模型对感染人数进行预测,结果如下图(下图是直接使用拟合得到的r值,为0.303,经过简单的测试,这个模型预测的数量偏小,如果有大佬对此比较了解能否指点一二):
2月9日预测结果(看最近几天的预测结果已经很接近了):
2月8日预测结果:
2月7日预测结果:
我们在高中生物中有学过,种群的增长曲线有分"J"型和"S"型,其中"J"型为理想型,实际中种群的数量不可能一直像指数函数那样,而是一个S型曲线,如下图所示。
这个函数的曲线和机器学习中的Sigmoid函数及其相似,其公式都有点像,logistic增长函数为”
其中K为环境最大容量,P0为初始容量,r为增长速率,r越大则增长越快(即更快的逼近上限)。
该模型的微分式是:
如果看不懂上面的内容,直接复制粘贴代码运行一下吧~(其实和之前的思路一样,使用最小二乘法拟合数据,只是拟合的函数换了而已)
from scipy.optimize import curve_fit
import urllib
import json
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import time
def date_encode(date):
# '01.24' -> 1 * 100 + 24 = 124
d = date.split('/')
month, day = int(d[0]), int(d[1])
return 100 * month + day
def date_decode(date):
# 124 -> '01.24'
return '{}.{}'.format(str(date // 100), str(date % 100))
def sequence_analyse(data):
date_list, confirm_list, dead_list, heal_list, suspect_list = [], [], [], [], []
data.sort(key = lambda x: date_encode(x['date']))
for day in data:
date_list.append(day['date'])
confirm_list.append(int(day['confirm']))
dead_list.append(int(day['dead']))
heal_list.append(int(day['heal']))
suspect_list.append(int(day['suspect']))
return pd.DataFrame({
'date': date_list,
'confirm': confirm_list,
'dead': dead_list,
'heal': heal_list,
'suspect': suspect_list
})
def get_date_list(target_month = 3):
"""
得到从1月13日到month月最后一天的所有日期列表
"""
month_day = [0, 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
cur_month, cur_day = 1, 13
ans = []
while cur_month <= target_month:
while cur_day <= month_day[cur_month]:
d = "0" + str(cur_day) if cur_day < 10 else str(cur_day)
ans += [str(cur_month) + "/" + d]
cur_day += 1
cur_day = 1
cur_month += 1
return ans
def logistic_function(t, K, P0, r):
t0 = 0
exp = np.exp(r * (t - t0))
return (K * exp * P0) / (K + (exp - 1) * P0)
def predict():
# 预测未来天数
predict_days = 20
# 获取实时数据
url = 'https://view.inews.qq.com/g2/getOnsInfo?name=wuwei_ww_cn_day_counts'
response = urllib.request.urlopen(url)
json_data = response.read().decode('utf-8').replace('\n','')
data = json.loads(json_data)
data = json.loads(data['data'])
df = sequence_analyse(data)
date, confirm = df['date'].values, df['confirm'].values
x = np.arange(len(confirm))
date_labels = get_date_list(4)
# 用最小二乘法估计拟合
popt, pcov = curve_fit(logistic_function, x, confirm)
print(popt)
#近期情况预测
predict_x = list(x) + [x[-1] + i for i in range(1, 1 + predict_days)]
predict_x = np.array(predict_x)
predict_y = logistic_function(predict_x, popt[0], popt[1], popt[2])
#绘图
plt.figure(figsize=(15, 8))
plt.plot(x, confirm, 's',label="confimed infected number")
plt.plot(predict_x, predict_y, 's',label="predicted infected number")
plt.xticks(predict_x, date_labels[:len(predict_x) + 1], rotation=90)
plt.yticks(rotation=90)
plt.suptitle("Logistic Fitting Curve for 2019-nCov infected numbers(Max = {}, r={:.3})".format(int(popt[0]), popt[2]), fontsize=16, fontweight="bold")
plt.title("Predict time:{}".format(time.strftime("%Y-%m-%d", time.localtime())), fontsize=16)
plt.xlabel('date', fontsize=14)
plt.ylabel('infected number', fontsize=14)
plt.plot()
predict()
参考: