前言

本文仅做根据已知点求拟合曲线的几种方法的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()

根据已知点求多项式

(以下内容来源网络）

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天感染人数:

总结

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
• 模型没有任何实践意义，且控制水平无法度量，预测无任何意义