寒暑假的梦幻联动可能是要过去了，最近的数模竞赛正好有关于新冠疫情的题目，于是想借机亲自分析这场梦幻联动的主角 COVID-19。
参数估计方法？
题目要求我们结合治愈率，死亡率等因素分析预测除中国外的某一国家的未来疫情数据，鉴于对于‘自由’的热爱，我们首先选择了自由美利坚来当做预测的目标，想要确定参数，我们首先需要做的就是寻找数据，本来自己想用python爬虫来解决，奈何自己爬虫的水平还是不够，最终还是在csdn上找了一篇代码想自己复制粘贴运行试试，没想到还真可以，万分感谢！后来看看代码发现其实也不难，主要是那个api网站的数据真是太便利了，下面就是那篇文章
疫情数据爬取
因此我们得到了美国共199天的疫情数据，接着我们使用这些数据来进行SEIR模型参数的估计，
传统的方程这里就不列了，这几个方程都烂大街了，，主要参数有三个，也就是易感者感染为潜伏者的概率、潜伏者变为感染者的概率，和移出感染者的概率，接着我们使用最小二乘法利用迭代方程，希望找出均方差最小的参数数值来得到这三个参数，我们在matlab中使用四层for循环（三个参数和199天的数据）来进行寻找参数的估计，我们想先用小范围，小步长的数据来进行逼近来防止程序的运行时间过长，后来发现循环后得到的均方差都很大，能达到10的14次方的数量级，因此我们寄希望于直接用小步长大范围来进行最优参数的寻找并进行了代码上的优化来减少循环的运行时间，不过那也运行了将近一个小时的时间，而且最后得到的误差也并不理想。
下面是matlab代码：
%--------------------------------------------------------------------------
% 初始化
x=1:199;%1月28-8月13共199天
t=1;
for i = 2:199
B(i,7)=(B(i,4)-B(i-1,4))/(B(i-1,2)-B(i-1,3)-B(i-1,4));
end
for i = 2:199
B(i,6)=(B(i,3)-B(i-1,3))/(B(i-1,2)-B(i-1,3)-B(i-1,4));
end
cure_rate=B(:,7)';%199天每天的治愈率数据
a1=polyfit(x,cure_rate,1);%得到拟合多项式函数
cure_rate1=polyval(a1,x);
dead_rate=B(:,6)';
a2=polyfit(x,dead_rate,1);%得到拟合多项式函数
dead_rate1=polyval(a2,x);
figure(1);
plot(x,cure_rate,'*',x,cure_rate1,'--r');
figure(2);
plot(x,dead_rate,'*',x,dead_rate1,'--b');
%--------------------------------------------------------------------------
% 参数设置
%--------------------------------------------------------------------------
N = 330000000; %人口总数
E = 0; %潜伏者
I = 5; %传染者
S = N - I; %易感者
R = 0; %康复者
D = 0; %死亡者
r=0.2; %感染率
%r = 1; %感染者接触易感者的人数
%b = 0.0003; %传染概率
r1 = 0.3; %潜伏者转化为感染者概率
r2=0.2;%潜伏者传染正常人概率
%r2 = 5; %潜伏者接触易感者的人数
%b2 = 0.03; %潜伏者传染正常人的概率
y = 0.01; %康复概率
d = 3e-4; %死亡概率
T = 1:199;
for r= 0:0.001:0.4
for r1 = 0:0.001:0.4
for r2 = 0:0.001:0.4
for idx = 1:199
% if polyval(a2,idx)<3e-5
% d=3e-5;
% else
% d=polyval(a2,idx);
% end
S(idx+1) = S(idx) - r*S(idx)*I(idx)/N(1) - r2*S(idx)*E(idx)/N;
E(idx+1) = E(idx) + r*S(idx)*I(idx)/N(1)-r1*E(idx) + r2*S(idx)*E(idx)/N;
I(idx+1) = I(idx) + r1*E(idx) - y*I(idx)-d*I(idx);
R(idx+1) = R(idx) + y*I(idx);
D(idx+1) = D(idx) + d*I(idx);
flag(idx) = (I(idx)-B(idx,2)-B(idx,3)-B(idx,4))^2;
end
error(t) = sum(flag);
t=t+1;
end
end
end
[v,address]=min(error);
address
v
在这之后我们分析了原因，误差很大的原因可能主要有三点，首先美国人口基数大，在均方误差上的数量级也更大，然后更为重要的原因是这疫情发生的199天中，各个国家的防控策略都在随着时间发生变化，而体现在美国的这199天的数据中，单一的参数并不适用，因此想要得到最为适合的模型参数，就需要对美国的这199天进行分阶段或者是用参数时变的思想进行处理，最后一点就是模型的问题，我们使用的是传统的SEIR模型，从这个模型方程中很容易发现，易感者-潜伏者-感染者-移出者这四个人群在进行单方向的流动，最后人群大概率都会变为移出者，因此这个传统的模型在国家层面上可能并不适用尤其是在把四个人群的总和确认为国家总人口时。
合理的参数估计
鉴于以上三点，我们放弃了研究美国的疫情走向（‘’自由’的锅？），而选择了目前疫情已经基本控制住的韩国，数据也是从上面的爬虫程序中得到的。考虑到上面的第二点，我们分析了我们的目标—即对某个国家未来的疫情走向做出分析，而并不需要考虑之前因防控策略不同所造成的影响，并观察韩国2月到8月的疫情走向图，因此我们果断选取了韩国近50天的数据（其趋势已经基本稳定下来），即下图中从150天到199天的数据。
然后考虑到上面的第三点，我们结合相关论文和资料认为这一单方向流动并不切合实际，并会给预测模型带来很大误差，因此我们在传统的SEIR模型中引入了潜伏者转为阴性的概率，并按照题目要求把移出者分为治愈者和死亡的人
之后也很自然的得到了迭代方程：
S(t+1) = S(t)-rbS(t)I(t)/N-r2b2S(t)E(t)/N+aE(t) E(t+1)=E(t)+rbS(t)I(t)/N-r1E(t)+r2b2S(t)E(t)/N-aE(t)
I(t+1)=I(t)+r1E(t)-yI(t)-dI(t);
R(t+1)=R(t)+yI(t)
D(t+1)=D(t)+dI(t)
下面是改进后SEIR模型中人群流动图
至于具体的参数估计我们考虑到四阶的循环有可能会得到较小误差，但是很可能并不符合实际，所以，我们参考下面的链接
日韩疫情参数估计结合现有形势得到感染率，潜伏者移出率等相关参数，另外，治愈率和死亡率，我们通过之前50天的治愈人数和死亡人数拟合得到。
最后我们得出模型曲线和实际数据曲线作对比，得到下图
可以从上图中看到模型估计感染者和实际感染者，估计康复者和实际康复者，估计死亡者和实际死亡者数据和走向上都吻合的很好，因此证明该模型及其参数具有一定的可信度。（当天晚上激动坏了QAQ）
下面是该模型的matlab代码
%--------------------------------------------------------------------------
% 初始化
load k.mat
x=1:50;%6月25-8月13共50天
t=1;
c=0;
for i = 2:50
k(i,7)=(k(i,3)-k(i-1,3))/(k(i-1,1)-k(i-1,2)-k(i-1,3));
end
for i = 2:50
k(i,6)=(k(i,2)-k(i-1,2))/(k(i-1,1)-k(i-1,2)-k(i-1,3));
end
cure_rate=k(:,7)';%199天每天的治愈率数据
a1=polyfit(x,cure_rate,1);%得到拟合多项式函数
cure_rate1=polyval(a1,x);
dead_rate=k(:,6)';
a2=polyfit(x,dead_rate,1);%得到拟合多项式函数
dead_rate1=polyval(a2,x);
figure(1);
plot(x,cure_rate,'*',x,cure_rate1,'--r');
figure(2);
plot(x,dead_rate,'*',x,dead_rate1,'--k');
%--------------------------------------------------------------------------
% 参数设置
%--------------------------------------------------------------------------
N = 51164435; %人口总数
E = 0; %潜伏者
I = 1307; %传染者
S = N - I; %易感者
R = 10974; %康复者
D = 282; %死亡者
r = 4.1176; %感染者接触易感者的人数
b = 0.017; %传染概率
r1 = 0.33; %潜伏者转化为感染者概率
r2 = 7; %潜伏者接触易感者的人数
b2 = 0.05; %潜伏者传染正常人的概率
y = 0.06; %康复概率
d = 0.5e-3; %死亡概率
T = 1:250;
for idx = 1:length(T)-1
if idx>=144
b=b*0.2;
b2=b2*0.2;
end
S(idx+1) = S(idx) - r*b*S(idx)*I(idx)/N(1) - r2*b2*S(idx)*E(idx)/N+0.5*E(idx);
E(idx+1) = E(idx) + r*b*S(idx)*I(idx)/N(1)-r1*E(idx) + r2*b2*S(idx)*E(idx)/N-0.5*E(idx);
I(idx+1) = I(idx) + r1*E(idx) - y*I(idx)-d*I(idx);
R(idx+1) = R(idx) + y*I(idx);
D(idx+1) = D(idx) + d*I(idx);
end
figure(3);
plot(T,E,T,I,T,R,T,D,1:50,k(:,1)-k(:,2)-k(:,3),'g*',1:50,k(:,2),'b*',1:50,k(:,3),'r*')
% plot(T,I,1:50,k(:,1)-k(:,2)-k(:,3))
% plot(1:50,k(:,1)-k(:,2)-k(:,3))
grid on;
xlabel('天');ylabel('人数')
% axis([0,70,0,5000])
legend('潜伏者','感染者','康复者','死亡者','实际感染人数','实际死亡人数','实际治愈人数')
%legend('感染者','实际感染人数');
%legend('实际感染人数');
for z= 68:7:250
c=c+1;
fn(c)=I(z);
end
figure(4)
plot(1:27,fn)
grid on;
xlabel('周');ylabel('感染者人数');
完。（疫情快过去吧！）
扫一扫
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
