前言
本文仅做根据已知点求拟合曲线的几种方法的python实现,无任何实际意义
数据来源(另一篇博文) 利用Python爬取新冠肺炎疫情实时数据,Pyecharts画2019-nCoV疫情地图
参考 https://blog.csdn.net/changdejie/article/details/83089933
数据读取
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
df = pd.read_excel("daily_data.xlsx")
df.head()
confirm | date | dead | heal | suspect | |
---|---|---|---|---|---|
0 | 41 | 1.13 | 1 | 0 | 0 |
1 | 41 | 1.14 | 1 | 0 | 0 |
2 | 41 | 1.15 | 2 | 5 | 0 |
3 | 45 | 1.16 | 2 | 8 | 0 |
4 | 62 | 1.17 | 2 | 12 | 0 |
根据已知点求多项式
(以下内容来源网络)
polyfit函数基于最小二乘法,使用的基本格式为:
p = polyfit(x,y,n)
[p,S] = polyfit(x,y,n)
[p,S,mu] = polyfit(x,y,n)
其中每个命令中的n为多项式拟合的次数,当n为1时,即为一次拟合(很多情况下等价于一元线性回归)。p是n+1维参数向量p(1),p(2)….那么拟合后对应的多项式即为p(1)*x^n + p(2)*x^(n-1) +…+ p(n)*x + p(n+1)。S是规模为1×1的结构数组,包括R(系数矩阵的QR分解的上三角阵),df(自由度),normr(拟合误差平方和的算术平方根)。
求出p之后我们需要作出拟合函数,那么只需要使用命令:
f=polyval(p,x)
然后plot出x和f即可。另外需要强调一点的是,往往需要在回归分析的时候给出相关系数(correlation coefficient),实际上也很简单,我们可以使用命令:
r=corrcoef(x,y);
这样得到的r即为相关系数矩阵,其中r(1,2)=r(2,1)为相关系数,其值在[-1,1]之间,1表示最大程度的正相关,-1表示最大程度的负相关。相关系数绝对值越靠近1,线性相关性质越好,根据数据描点画出来的函数-自变量图线越趋近于一条平直线,拟合的直线与描点所得图线也更相近。
另外,转载两条使用polyfit的注意事项:
1.使用polyfit命令进行多项式拟合时要注意的是,向量x(其中元素作为自变量)中不重复的元素个数m,和拟合阶数k需要满足m>=k+1.简单分析:k阶拟合需要确定k+1个未知参数(如1阶拟合y = ax + b需要确定a和b两个参数),故而至少需要k+1个方程,故而需要至少k+1个不同的已知数对(x,y),由于函数中x只能对应一个y,故而需要至少k+1个不同的x。
2.polyfit只适合于形如y = a[k]*x^k + a[k-1]*x^(k-1) + …. + a[1]*x + a[0]的完全的一元多项式的数据拟合。
t = list(range(13,13 + len(df),1))
t = np.array(t).astype(int)
y = np.array(list(df["confirm"])).astype(int)
print("时间t :\n",t)
print("确诊人数 :\n",y)
f1 = np.polyfit(t,y,3)
print("a,b,c,d = :\n",f1)
# 返回拟合函数
p1 = np.poly1d(f1)
print("多项式曲线拟合函数 :\n",p1)
yvals = np.polyval(f1,t)
plot1 = plt.plot(t,y,"ob")
plot2 = plt.plot(t,yvals,"r")
plt.xlabel("t (t0=13)")
plt.ylabel("P_confirm")
plt.show()
时间t :
[13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29]
确诊人数 :
[ 41 41 41 45 62 198 275 291 440 571 830 1287 1975 2744
4515 5974 7711]
a,b,c,d = :
[ 5.26169591e+00 -2.76422472e+02 4.80349024e+03 -2.74304487e+04]
多项式曲线拟合函数 :
3 2
5.262 x - 276.4 x + 4803 x - 2.743e+04
马儿萨斯模型-指数增长(J型)
模型解释
-
问题1:在知道当前或过去某个时刻的人口数量的情况下,预测未来某个时刻的人口数量?
-
问题2:t趋于无穷时,即在遥远的未来,人口趋势是怎样的?
模型求解
#定义函数
def L(t,p0,r1):
return p0*np.exp(r1*t)
popt,pcov = curve_fit(L,t,y)
P0 = popt[0]
r1 = popt[1]
yvals = L(t,P0,r1)
print("按照J型增长,自然增长率 r1 = " + str(r1) + " (天)推算,从发病之日起(12月31日):")
plot1 = plt.plot(t,y,"p")
plot2 = plt.plot(t,yvals,"m")
for i in range(30,80,5):
if int(L(i,P0,r1)) <100000000000:
print(str(i) + "天后,预计感染" + str(int(L(i,P0,r1))))
按照J型增长,自然增长率 r1 = 0.34243812944795793 (天)推算,从发病之日起(12月31日):
30天后,预计感染11286
35天后,预计感染62538
40天后,预计感染346529
45天后,预计感染1920151
50天后,预计感染10639727
55天后,预计感染58955655
60天后,预计感染326678400
65天后,预计感染1810153360
70天后,预计感染10030216814
75天后,预计感染55578301571
阻滞增长模型(S型)
模型解释
- X(t) t时刻的感染数量
- X0 初始感染数量
- r 每日净增长率
- Xm 环境所能容纳的最大感染数量
模型求解
def S(t,K,p0,r2):
t0 = 13 #定义13号为初始时间
exp_value=np.exp(r2*(t-t0))
return (K*exp_value*p0)/(K+(exp_value-1)*p0)
popt,pcov = curve_fit(S,t,y)
K = popt[0]
p0 = popt[1]
r2 = popt[2]
print("按照目前的控制水平,阻滞增长模型S型增长率为" + str(r2) + ",最大患病数为" +str(int(K)))
yvals_1 = S(t,K,p0,r2)
plot1 = plt.plot(t,y,"ob")
plot1 = plt.plot(t,yvals_1,"r")
按照目前的控制水平,阻滞增长模型S型增长率为0.47179446312327517,最大患病数为15870
print("按照29号之前的控制方案,未来10天感染人数:")
data = []
for i in range(30,40):
if i < 31:
data.append(["1-" + str(i + 1),int(S(i,K,p0,r2))])
else:
data.append(["2-" + str(i-int(30)),int(S(i,K,p0,r2))])
pre_data = pd.DataFrame(data)
pre_data
按照29号之前的控制方案,未来10天感染人数:
0 | 1 | |
---|---|---|
0 | 1-31 | 9606 |
1 | 2-1 | 11281 |
2 | 2-2 | 12657 |
3 | 2-3 | 13700 |
4 | 2-4 | 14443 |
5 | 2-5 | 14949 |
6 | 2-6 | 15282 |
7 | 2-7 | 15498 |
8 | 2-8 | 15636 |
9 | 2-9 | 15723 |
总结
t1 = np.array(range(13,60,1))
model_F = np.polyval(f1,t1)
model_L = L(t1,P0,r1)
model_S = S(t1,K,p0,r2)
plt.figure(figsize=(12, 8))
plt.subplot(2, 2, 1)
plot1 = plt.plot(t,y,"b")
plt.title('Actual diagnosis')
plt.subplot(2, 2, 2)
plot2 = plt.plot(t1,model_F,"m")
plt.title('Polynomial prediction')
plt.subplot(2, 2, 3)
plot3 = plt.plot(t1,model_L,"k")
plt.title('Model_L prediction')
plt.subplot(2, 2, 4)
plot4 = plt.plot(t1,model_S,"r")
plt.title('Model_S prediction')
几点总结(无实际意义)
- 无任何控制措施的情况下(J型),预计从12-31日起,65天-70天全球感染
- 按照29号之前的控制措施(S型),预计2-15日左右达到最大感染量15870
- 模型没有任何实践意义,且控制水平无法度量,预测无任何意义