年初的新冠疫情来势汹汹,但好在政府及时控制住,经济得以恢复正常。疫情发生后,国内外很多研究学者都通过建模等方法分析了疫情可能导致的感染人数,下面分享一下通过Matlab的多项式曲线拟合预测新冠病毒感染人数趋势,结果粗糙,仅仅作为学习。
写在前言:预测数据基于官方公布的数据,本文只做一个学术研究,结果仅供参考;
1、数据获取
为了简化过程,网上查询到对应的感染人数。选择了2020/1/18到2/3日的数据
x1 = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]; % 累计天数
x2 = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]; % 新增天数
y1 = [291,391,440,571,830,1287,1975,2744,4515,5974,7711,9692,11791,14411,17205,20438]; % 累计
y2 = [77,149,131,259,444,688,769,1771,1459,1737,1982,2102,2590,2829,3235]; % 新增
2、数据处理
在后续绘图的过程中,将横坐标改为日期的形式,用来表示每天的病例人数,用到了matlab中的日历类型。
首先,使用datetime函数,指定开始日期:
t1=datetime(2020,1,18);
由于绘图时横坐标需要精确到每天,故将每天的datetime类型通过for循环存入数组:
利用caldays函数进行日期的加减运算。
t1=datetime(2020,1,18);
for i = 1:length(x1)
accumulate = [accumulate t1+caldays(i)];
end
for i = 1:length(x2)
add = [add t1+caldays(i+1)];
end
for i=1:33
date = [date t1+caldays(i)];
end
这样,使用第二段代码中data这个数组作为横坐标就能绘制日期横坐标的图了。
3、模型建立——多项式拟合
3.1、多项式拟合原理和本文说明
多项式拟合是用一个多项式展开去拟合包含数个分析格点的一小块分析区域中的所有观测点,得到观测数据的客观分析场。展开系数用最小二乘拟合确定。
多项式拟合,将一组数据尽可能的映射到一个多项式函数上,反映这一组数据之间的一个函数关系。故可以使用多项式拟合方法对于数据的未来走向进行一定程度上的预测。
对多项式进行曲线拟合可以使用polyfit函数,该函数能够很好地进行曲线拟合,用法MATLAB程序代码为:
p =polyfit(x,y,n)
其中,x为横坐标,在本文中,为自2020-01-18开始的天数(1代表2020-01-19,以此类推),y为纵坐标,在本为中,为自2020-01-18开始每日的累计感染人数。
n为拟合的多项式次数。
3.2、拟合次数n的选取
在本文中,n=4,即多项式的最高次数为4。选取原因如下:
1、一次、二次、三次函数单调性太过明显,线性和指数级别的单调增加不能良好反映疫情的发展状况;
2、五次以上的拟合情况,结果如下:
据分析,n=3时,单调递增显然不符合疫情发展趋势,n=5和6时结果太夸张。
所以3、5、6和更高次的拟合都不符合目前的实际情况。
当n=4时:
在四次多项式拟合中疫情出现了拐点,虽然在拐点之后,函数是单调递减的,但是结合目前国家的大力防控措施,可以姑且认为四次多项式在拐点之前的增长是具有一定参考价值的。
拟合代码实现如下(使用polyval函数预测每日感染人数):
order = 4;
p1 =polyfit(x1, y1, order); % 累计 四次函数
p2 =polyfit(x2, y2, order-1); % 新增是累计的一阶导数 三次函数
add_new = polyval(p2,day_num); % 预测的新增
total = polyval(p1,day_num); % 累计
画图函数实现如下:
figure
plot(date,total,'r')
hold on;
plot(date,add_new,'b')
hold on;
plot(accumulate, y1, '*')
hold on;
plot(add, y2, 'o')
grid on;
4、数据分析与预测
通过将1月18日统计数据开始后的天数代入函数,使用polyval函数计算每日的感染人数,并绘制图像
由此分析,疫情将在2月15日至2月20日期间出现“拐点”,即累计感染人数显著降低。
代码全文:
clc
clear
close all
x1 = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]; % 累计天数
x2 = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]; % 新增天数
y1 = [291,391,440,571,830,1287,1975,2744,4515,5974,7711,9692,11791,14411,17205,20438]; % 累计
y2 = [77,149,131,259,444,688,769,1771,1459,1737,1982,2102,2590,2829,3235]; % 新增
figure
plot(x1,y1,x2,y2)
xlabel 天
ylabel 人数
legend('累计人数','新增人数')
accumulate = []; % 累计和新增天数横坐标
add = [];
date = []; % 累计和新增天数横坐标
day_num = 1:33; % 预测的天数;
t1=datetime(2020,1,18);
for i = 1:length(x1)
accumulate = [accumulate t1+caldays(i)];
end
for i = 1:length(x2)
add = [add t1+caldays(i+1)];
end
for i=1:33
date = [date t1+caldays(i)];
end
figure
plot(accumulate,y1,add,y2)
xlabel 天
ylabel 人数
legend('累计人数','新增人数')
order = 4;
p1 =polyfit(x1, y1, order); % 累计 四次函数
p2 =polyfit(x2, y2, order-1); % 新增是累计的一阶导数 三次函数
add_new = polyval(p2,day_num); % 预测的新增
total = polyval(p1,day_num); % 累计
figure
plot(date,total,'r')
hold on;
plot(date,add_new,'b')
hold on;
plot(accumulate, y1, '*')
hold on;
plot(add, y2, 'o')
grid on;
% title('拟合次数为4')
最后祝各位读者朋友2020远离疾病、远离灾难,平安渡过2020