本文首先介紹了Logistic模型的原理，然后嘗試用Logistic曲線擬合疫情攒暇，雖然疫情已經(jīng)接近尾聲土匀，模型的預(yù)測意義并不大，但仍然可以通過回溯過去發(fā)現(xiàn)有趣的現(xiàn)象形用。
目錄
1. Logistic模型
??1.1 馬爾薩斯人口模型
??1.2 Logistic增長模型
2. 用Logistic模型擬合疫情
??2.1 湖北省
??2.2 湖北省外
3. Python實(shí)現(xiàn)
參考文獻(xiàn)
1. Logistic模型
邏輯斯蒂方程(Logistic function)由比利時(shí)數(shù)學(xué)家兼生物學(xué)家皮埃爾·弗朗索瓦·韋呂勒(Pierre Francois Verhulst)在研究人口增長模型時(shí)提出就轧，是對馬爾薩斯人口模型(Malthus, 1798)的改進(jìn)。下圖可以直觀看出兩者的區(qū)別田度，下面將對兩種模型作詳細(xì)介紹妒御。
1.1 馬爾薩斯人口模型
馬爾薩斯人口模型假定人口增長率保持不變，
其中表示人口數(shù)镇饺，是時(shí)間
的函數(shù)乎莉。
求解微分方程可以得到人口隨時(shí)間變化的函數(shù)，
其中表示第0期人口。
不難發(fā)現(xiàn)惋啃，人口呈現(xiàn)指數(shù)增長哼鬓，即J型曲線。然而現(xiàn)實(shí)中受到自然資源約束以及疾病等因素影響边灭，人口增長率不可能一直保持不變魄宏。
1.2 Logistic增長模型
皮埃爾在馬爾薩斯人口模型的基礎(chǔ)上進(jìn)行了改進(jìn)，將人口增長率設(shè)為存筏，其中
可以理解為環(huán)境最大允許的最大人口數(shù)量宠互。此時(shí)，當(dāng)人口
越接近于
時(shí)椭坚，增長率越低予跌，即人口增長率隨人口數(shù)量的增加而線性減少。
求解微分方程可以得到人口隨時(shí)間變化的函數(shù)善茎，
其中表示第0期人口券册。
如下左圖，該曲線描述的人口增長呈現(xiàn)S型垂涯，增長速率隨時(shí)間先增后減烁焙，在處增長最快。注意到耕赘，增長速率(
)表示人口當(dāng)期變化的絕對數(shù)值骄蝇，而增長率(
)表示人口變化量與當(dāng)期人口的比值。
相比于馬爾薩斯人口模型操骡，Logistic增長模型更加符合實(shí)際九火，該模型常常被應(yīng)用于描述種群、傳染病增長以及商品銷售量預(yù)測等領(lǐng)域册招。
2. 用Logistic模型擬合疫情
考慮到湖北省內(nèi)與湖北省外存在異質(zhì)性岔激，下面將用Logistic模型分別對湖北省與湖北省外的累計(jì)確診人數(shù)進(jìn)行擬合。首先需要獲取累計(jì)確診人數(shù)是掰，數(shù)據(jù)來源為WindQuant提供的Wind數(shù)據(jù)接口虑鼎，更多數(shù)據(jù)下載方法見獲取COVID-19疫情歷史數(shù)據(jù)的n種方法。
然后使用Scipy.optimezi庫的curve_fit函數(shù)對Logistic曲線進(jìn)行非線性最小二乘擬合键痛。待定參數(shù)包括三個(gè)散休，根據(jù)最小化MSE準(zhǔn)則媒楼，采用網(wǎng)格調(diào)參的方式尋找最優(yōu)參數(shù)：對最大容量
以步長1遍歷（10000, 80000）區(qū)間，對增長率
以步長0.01遍歷（0, 1）區(qū)間戚丸。
2.1 湖北省
湖北省疫情擬合結(jié)果：最大容量為67667划址，增長率
為0.24扔嵌。疫情的拐點(diǎn)發(fā)生在2月8日前后，這也是疫情增長最快一段時(shí)期夺颤，在這段時(shí)間內(nèi)痢缎，由于疫情的快速增長和醫(yī)療資源的相對匱乏，官方公布的確診人數(shù)失真嚴(yán)重世澜。從擬合結(jié)果來看独旷，當(dāng)前疫情已經(jīng)到達(dá)尾聲。
2.2 湖北省外
省外疫情擬合結(jié)果：最大容量為12857寥裂，增長率
為0.26嵌洼。疫情的拐點(diǎn)發(fā)生在2月3日前后，當(dāng)前疫情也已經(jīng)到達(dá)尾聲封恰。
3. Python實(shí)現(xiàn)
下面直接給出Python代碼麻养，注意第17行需要將userid修改為自己的WindQuant用戶ID，當(dāng)然也可以對代碼略作修改后讀取本地?cái)?shù)據(jù)或者其他數(shù)據(jù)源诺舔。微信搜索公眾號PurePlay鳖昌，后臺回復(fù)Logistic，即可獲取完整的數(shù)據(jù)低飒、代碼以及結(jié)果文件许昨。
import numpy as np
import pandas as pd
import datetime as dt
import time
import requests
import json
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams['font.sans-serif'] = ['SimHei'] # 用來正常顯示中文標(biāo)簽
mpl.rcParams['axes.unicode_minus'] = False # 用來正常顯示負(fù)號
from scipy.optimize import curve_fit
from sklearn.metrics import mean_squared_error
def load_data():
# 下載原始數(shù)據(jù)，數(shù)據(jù)來源為WindQuant提供的Wind數(shù)據(jù)接口
# "Nation","Hubei","Outside"褥赊，分別表示全國糕档、省內(nèi)和省外的累計(jì)確診人數(shù)
userid = "94ac606a-c220-46af-abf3" # 此處僅為無效id，需要改為自己的WindQuant用戶ID
indicators = "S6274770,S6275447,S6291523"
factors_name = ["Nation","Hubei","Outside"]
startdate = "2020-01-16"
enddate = dt.datetime.strftime(dt.date.today() + dt.timedelta(-1),"%Y-%m-%d")# 數(shù)據(jù)的結(jié)束日期設(shè)置為昨天
url = '''https://www.windquant.com/qntcloud/data/edb?userid={}&indicators={}&startdate={}&enddate={}'''.format(
userid,indicators,startdate,enddate)
response = requests.get(url)
data = json.loads(response.content.decode("utf-8"))
try:
time_list = data["times"]
value_list = data["data"]
for i in range(len(time_list)):
time_list[i] = time.strftime("%Y-%m-%d", time.localtime(time_list[i]/1000))
result = pd.DataFrame(columns=factors_name, index = time_list)
for i in range(len(factors_name)):
result[factors_name[i]] = value_list[i]
print(result)
result.to_csv(r"Data\LogisticData.csv")
return result
except Exception as e:
print("服務(wù)異常")
def data_abstract(result,area):
# result：下載的原始數(shù)據(jù)
# area：“Hubei”或“Outside”崭倘，分別返回省內(nèi)和省外的累計(jì)確診人數(shù)
y_data = result[area]
# 刪除缺失值并轉(zhuǎn)換為float型
y_data[y_data == 'NaN'] = np.NAN
y_data = y_data.dropna()
y_data = y_data.astype(float)
global first_date # 后續(xù)數(shù)據(jù)可視化需要
first_date = dt.datetime.strptime(y_data.index[0],'%Y-%m-%d')
# 返回與數(shù)據(jù)集等長的從0開始的時(shí)間序列作為logistic函數(shù)的自變量
x_data = np.asarray(range(0,len(y_data)))
return x_data, y_data
hyperparameters_r = None
hyperparameters_K = None
def logistic_increase_function(t,P0):
# logistic生長函數(shù)：t:time P0:initial_value K:capacity r:increase_rate
# 后面將對r和K進(jìn)行網(wǎng)格優(yōu)化
r = hyperparameters_r
K = hyperparameters_K
exp_value = np.exp(r * (t))
return (K * exp_value * P0) / (K + (exp_value - 1) * P0)
def fitting(logistic_increase_function, x_data, y_data):
# 傳入要擬合的logistic函數(shù)以及數(shù)據(jù)集
# 返回?cái)M合結(jié)果
popt = None
mse = float("inf")
i = 0
# 網(wǎng)格搜索來優(yōu)化r和K參數(shù)
r = None
k = None
k_range = np.arange(10000, 80000, 1000)
r_range = np.arange(0, 1, 0.01)
for k_ in k_range:
global hyperparameters_K
hyperparameters_K = k_
for r_ in r_range:
global hyperparameters_r
hyperparameters_r = r_
# 用非線性最小二乘法擬合
popt_, pcov_ = curve_fit(logistic_increase_function, x_data, y_data, maxfev = 4000)
# 采用均方誤準(zhǔn)則選擇最優(yōu)參數(shù)
mse_ = mean_squared_error(y_data, logistic_increase_function(x_data, *popt_))
if mse_ <= mse:
mse = mse_
popt = popt_
r = r_
k = k_
i = i+1
print('\r當(dāng)前進(jìn)度：{0}{1}%'.format('▉'*int(i*10/len(k_range)/len(r_range)),int(i*100/len(k_range)/len(r_range))), end='')
print('擬合完成')
hyperparameters_K = k
hyperparameters_r = r
popt, pcov = curve_fit(logistic_increase_function, x_data, y_data)
print("K:capacity P0:initial_value r:increase_rate")
print(hyperparameters_K, popt, hyperparameters_r)
return hyperparameters_K, hyperparameters_r, popt
def predict(logistic_increase_function, popt):
# 根據(jù)最優(yōu)參數(shù)進(jìn)行預(yù)測
future = np.linspace(0, 60, 60)
future = np.array(future)
future_predict = logistic_increase_function(future, popt)
diff = np.diff(future_predict)
diff = np.insert(diff, 0, np.nan)
return future, future_predict, diff
def visualize(area, future, future_predict, x_data, y_data, diff):
# 繪圖
x_show_data_all = [(first_date + (dt.timedelta(days=i))).strftime("%m-%d") for i in future]
x_show_data = x_show_data_all[:len(x_data)]
plt.figure(figsize=(12, 6), dpi=300)
plt.scatter(x_show_data, y_data, s=35, marker='.', c = "dimgray", label="確診人數(shù)")
plt.plot(x_show_data_all, future_predict, 'r', linewidth=1.5, label='預(yù)測曲線')
plt.plot(x_show_data_all, diff, "r", c='darkorange',linewidth=1.5, label='一階差分')
plt.tick_params(labelsize=10)
plt.xticks(x_show_data_all)
plt.grid() # 顯示網(wǎng)格
plt.legend() # 指定legend的位置右下角
ax = plt.gca()
for label in ax.xaxis.get_ticklabels():
label.set_rotation(45)
if area == "Hubei":
plt.ylabel('湖北省累計(jì)確診人數(shù)')
plt.savefig(r"Data\Hubei.png",dpi=300)
else:
plt.ylabel('湖北省外累計(jì)確診人數(shù)')
plt.savefig(r"Data\Outside.png",dpi=300)
plt.show()
if __name__ == '__main__':
# 載入數(shù)據(jù)
result = load_data()
for area in ["Outside", "Hubei"]:
# 從原始數(shù)據(jù)中提取對應(yīng)數(shù)據(jù)
x_data, y_data = data_abstract(result, area)
# 擬合并通過網(wǎng)格調(diào)參尋找最優(yōu)參數(shù)
K, r, popt = fitting(logistic_increase_function, x_data, y_data)
# 模型預(yù)測
future, future_predict, diff= predict(logistic_increase_function, popt)
# 繪制圖像
visualize(area, future, future_predict, x_data, y_data, diff)
參考文獻(xiàn)
[1] Verhulst, P. E., Corresp. Math. Phys., 10, 113, 1838.
[2] 馬爾薩斯．人口原理[M]．北京：商務(wù)印書館翼岁，1992.6．
[2] 宋波, 玄玉仁, 盧鳳勇, et al. 淺評邏輯斯蒂方程[J]. 生態(tài)學(xué)雜志, 1986(03):59-64.
[3] 徐榮輝. 邏輯斯蒂方程及其應(yīng)用[J]. 山西財(cái)經(jīng)大學(xué)學(xué)報(bào),2010,32(S2):311-312.
相關(guān)文章
獲取COVID-19疫情歷史數(shù)據(jù)的n種方法
AHP | 層次分析法原理及Python實(shí)現(xiàn)
以上是本篇的全部內(nèi)容，歡迎關(guān)注我的知乎|簡書|CSDN|微信公眾號PurePlay , 會不定期分享研究與學(xué)習(xí)干貨司光。
