疫情数据分析和预测

疫情数据分析和预测是医学和流行病学应对大范围流行病时的重要判断手段，在医治隔离、预防响应、物资生产调配等抗疫措施上起到参考作用。

以下将通过已知模型尝试寻找合适拟合模型并对目前全球疫情发展作出一定程度的预测。

一、逻辑斯蒂模型（Logistic）

（1）模型描述：当一个物种迁入到一个新生态系统中后，其数量会发生变化。假设该物种的起始数量小于环境的最大容纳量，则数量会增长。该物种在此生态系统中有天敌、食物、空间等资源也不足（非理想环境），则增长函数满足逻辑斯谛方程，图像呈S形，此方程是描述在资源有限的条件下种群增长规律的一个最佳数学模型。

（2）一般疾病的传播是S型增长的过程，因为疾病传播的过程中会受到一定的阻力（医治、切断传播途径等措施）。

此处采用最小二乘法，对logistic增长函数进行拟合。以下将检验最小二乘法拟合的逻辑斯蒂模型是否能贴合实际。

import numpy as np
import matplotlib.pyplot as plt
import math
import pandas as pd
 
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
 
def logistic_increase_function(t,K,P0,r):
    t0=11
    #r 0.05/0.55/0.65
    r = 0.45
    # t:time   t0:initial time    P0:initial_value    K:capacity  r:increase_rate
    exp_value=np.exp(r*(t-t0))
    return (K*exp_value*P0)/(K+(exp_value-1)*P0)
 
'''
1.11日41例
1.18日45例
1.19日62例
1.20日291例
1.21日440例
1.22日571例
1.23日835例
1.24日1297例
1.25日1985例
1.26日2762例
1.27日4535例
'''
 
#  日期及感染人数
t=[11,18,19,20 ,21, 22, 23, 24,  25,  26,  27]
t=np.array(t)
P=[41,45,62,291,440,571,835,1297,1985,2762,4535]
P=np.array(P)
 
# 用最小二乘法估计拟合
# 现有数据曲线拟合检验
popt1, pcov1 = curve_fit(logistic_increase_function, t, P)

#获取popt里面是拟合系数
print("K:capacity  P0:initial_value   r:increase_rate   t:time")
print(popt1)
#拟合后预测的P值

P_predict = logistic_increase_function(t,popt1[0],popt1[1],popt1[2])
#未来预测
future=[11,18,19,20 ,21, 22, 23, 24,  25,  26,  27,28,29,30,31,41,51,61,71,81,91,101]
future=np.array(future)
future_predict=logistic_increase_function(future,popt1[0],popt1[1],popt1[2])

#近期情况预测
tomorrow=[28,29,30,32,33,35,37,40]
tomorrow=np.array(tomorrow)
tomorrow_predict=logistic_increase_function(tomorrow,popt1[0],popt1[1],popt1[2])
#绘图
plot1 = plt.plot(t, P, 's',label="confimed infected people number")
plot2 = plt.plot(t, P_predict, 'r',label='predict infected people number')
plot3 = plt.plot(tomorrow, tomorrow_predict, 's',label='tomorrow 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),popt1[0],popt1[1],popt1[2]))
print(logistic_increase_function(np.array(29),popt1[0],popt1[1],popt1[2]))
plt.show()

本次拟合采用了1月11日到1月27日的累计确诊病例数据作为原始数据，采用最小二乘法拟合逻辑斯蒂曲线，最后经过对逻辑斯蒂模型中R值（增长速率，到达K值的速度）的拟合调整，发现在0.45附近得到的曲线比较贴合我国1月至2月疫情实际情况。2月9日的预测值在4万左右，与实际情况十分贴近，也证明了模型的一定可靠性。

将本模型推广，进行全球范围内典型新冠肺炎爆发国家的疫情拟合与未来疫情预测，同时将通过R值的大小反应出该国疫情应对的有效程度。

import numpy as np
import matplotlib.pyplot as plt
import math
import pandas as pd
 
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
 
def logistic_increase_function(t,K,P0,r):
    t0=11
    #r 
    r = 0.08
    # t:time   t0:initial time    P0:initial_value    K:capacity  r:increase_rate
    exp_value=np.exp(r*(t-t0))
    return (K*exp_value*P0)/(K+(exp_value-1)*P0)
  
#  日期及感染人数

n = "../dataSets\\countrydata.csv"
data = pd.read_csv(n)
# 修改国家可以得到不同的曲线拟合情况
data = data[data['countryName'] == '日本']
date_list = list(data['dateId'])
date_list = list(map(lambda x:str(x),date_list))

confirm_list = list(data['confirmedCount'])

time_array = np.array(range(19,len(date_list)+19))
long_time_array = np.array(range(19,len(date_list)+190))
confirm_array = np.array(confirm_list)
 
# 用最小二乘法估计拟合

# 现有数据曲线拟合预测
popt, pcov = curve_fit(logistic_increase_function, time_array, confirm_array)
#获取popt里面是拟合系数
print("K:capacity  P0:initial_value   r:increase_rate   t:time")
print(popt)
#拟合后预测的P值
P_predict = logistic_increase_function(long_time_array,popt[0],popt[1],popt[2])

#未来预测

#近期情况预测

#绘图
plot1 = plt.plot(time_array, confirm_array, 's',label="confimed infected people number")
plot2 = plt.plot(long_time_array, P_predict, 'r',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()

import numpy as np
import matplotlib.pyplot as plt
import math
import pandas as pd
 
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
 
def logistic_increase_function(t,K,P0,r):
    t0=11
    #r 中国0.25 美国0.05 英国0.08 意大利0.08 德国0.09 韩国0.11
    r = 0.05
    # t:time   t0:initial time    P0:initial_value    K:capacity  r:increase_rate
    exp_value=np.exp(r*(t-t0))
    return (K*exp_value*P0)/(K+(exp_value-1)*P0)
 
#  日期及感染人数

n = "../dataSets\\countrydata.csv"
data = pd.read_csv(n)
# 修改国家可以得到不同的曲线拟合情况
data = data[data['countryName'] == '美国']
date_list = list(data['dateId'])
date_list = list(map(lambda x:str(x),date_list))

confirm_list = list(data['confirmedCount'])

time_array = np.array(range(19,len(date_list)+19))
long_time_array = np.array(range(19,len(date_list)+190))
confirm_array = np.array(confirm_list)
 
# 用最小二乘法估计拟合

# 现有数据曲线拟合预测
popt, pcov = curve_fit(logistic_increase_function, time_array, confirm_array)
#获取popt里面是拟合系数
print("K:capacity  P0:initial_value   r:increase_rate   t:time")
print(popt)
#拟合后预测的P值
P_predict = logistic_increase_function(long_time_array,popt[0],popt[1],popt[2])

#未来预测

#近期情况预测

#绘图
plot1 = plt.plot(time_array, confirm_array, 's',label="confimed infected people number")
plot2 = plt.plot(long_time_array, P_predict, 'r',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()

import numpy as np
import matplotlib.pyplot as plt
import math
import pandas as pd
 
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
 
def logistic_increase_function(t,K,P0,r):
    t0=11
    #r 
    r = 0.09
    # t:time   t0:initial time    P0:initial_value    K:capacity  r:increase_rate
    exp_value=np.exp(r*(t-t0))
    return (K*exp_value*P0)/(K+(exp_value-1)*P0)
 
#  日期及感染人数

n = "../dataSets\\countrydata.csv"
data = pd.read_csv(n)
# 修改国家可以得到不同的曲线拟合情况
data = data[data['countryName'] == '德国']
date_list = list(data['dateId'])
date_list = list(map(lambda x:str(x),date_list))

confirm_list = list(data['confirmedCount'])

time_array = np.array(range(19,len(date_list)+19))
long_time_array = np.array(range(19,len(date_list)+190))
confirm_array = np.array(confirm_list)
 
# 用最小二乘法估计拟合

# 现有数据曲线拟合预测
popt, pcov = curve_fit(logistic_increase_function, time_array, confirm_array)
#获取popt里面是拟合系数
print("K:capacity  P0:initial_value   r:increase_rate   t:time")
print(popt)
#拟合后预测的P值
P_predict = logistic_increase_function(long_time_array,popt[0],popt[1],popt[2])

#未来预测

#近期情况预测

#绘图
plot1 = plt.plot(time_array, confirm_array, 's',label="confimed infected people number")
plot2 = plt.plot(long_time_array, P_predict, 'r',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()

import numpy as np
import matplotlib.pyplot as plt
import math
import pandas as pd
 
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
 
def logistic_increase_function(t,K,P0,r):
    t0=11
    #r 中国0.25 美国0.05 英国0.08 意大利0.08 德国0.09 韩国0.11
    r = 0.11
    # t:time   t0:initial time    P0:initial_value    K:capacity  r:increase_rate
    exp_value=np.exp(r*(t-t0))
    return (K*exp_value*P0)/(K+(exp_value-1)*P0)
 
#  日期及感染人数

n = "../dataSets\\countrydata.csv"
data = pd.read_csv(n)
# 修改国家可以得到不同的曲线拟合情况
data = data[data['countryName'] == '韩国']
date_list = list(data['dateId'])
date_list = list(map(lambda x:str(x),date_list))

confirm_list = list(data['confirmedCount'])

time_array = np.array(range(19,len(date_list)+19))
long_time_array = np.array(range(19,len(date_list)+190))
confirm_array = np.array(confirm_list)
 
# 用最小二乘法估计拟合

# 现有数据曲线拟合预测
popt, pcov = curve_fit(logistic_increase_function, time_array, confirm_array)
#获取popt里面是拟合系数
print("K:capacity  P0:initial_value   r:increase_rate   t:time")
print(popt)
#拟合后预测的P值
P_predict = logistic_increase_function(long_time_array,popt[0],popt[1],popt[2])

#未来预测

#近期情况预测

#绘图
plot1 = plt.plot(time_array, confirm_array, 's',label="confimed infected people number")
plot2 = plt.plot(long_time_array, P_predict, 'r',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()