1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
| %% SIR差分模型
frame=importdata('美国covid19疫情数据l-history改.csv');
date_row=size(frame.data,1);
data=frame.data;
date=frame.textdata(:,1);
real_date=date(2:date_row+1,:);
E=zeros(2,90);
l=0.000125; %日接触率
m=0.01; %日治愈率0.0005
E(1,1)=data(345,12)/300000000;
E(2,1)=1-E(1,1);
for i=1:90
E(1,i+1)=l*E(2,i)-m*E(1,i)+E(1,i);
E(2,i+1)=E(2,i)-l*E(1,i)*E(2,i);
end
X=flip(data(311:345,12)');% 5.31 4.2
rate=flip(data(311:345,12)')%现阳性率
rate(isnan(rate))=0;
n=length(rate);
rt=0:1:n-1;
%a0=[100,10]; %初值
figure(1)
h1 = plot(rt,X/300000000,'*'); %画点
hold on;
a=E(1,:);
b=E(2,:);
h2 = plot(a,'r')
xlabel('日期'); %设置横坐标名
ylabel('感染人数'); %设置纵坐标名
legend([h1 h2],'3~5月感染人数','SIR模型迭代曲线','Location','NorthWest');
hold on
%% 预测
X=flip(data(255:345,12)');%311行后为第一题
rate=flip(data(255:345,12)')%现阳性率
rate(isnan(rate))=0;
n=length(rate);
rt=0:1:n-1;
figure(2)
h3 = plot(rt,X/300000000,'*'); %画点
hold on;
a=E(1,:);
b=E(2,:);
h4=plot(a,'r')
xlabel('日期');
ylabel('感染人数');
legend([h3 h4],'3~6月感染人数','SIR模型迭代曲线','Location','NorthWest');
hold on
%% 接种疫苗前
frame=importdata('usc.csv');
date_row=size(frame.data,1);
data=frame.data;
date=frame.textdata(:,1);
real_date=date(2:date_row+1,:);
E=zeros(3,150);
l=0.0018; %日接触率 0.0018
m=0.01; %日治愈率 0.01
E(1,1)=data(320,1)/300000000;
E(3,1)=28664448/300000000%43714928 33714928 28664448
E(2,1)=1-E(1,1)-E(3,1);
for i=1:150
E(1,i+1)=E(1,i)+l*E(2,i)-m*E(1,i);
E(2,i+1)=E(2,i)-l*E(1,i)*E(2,i);
E(3,i+1)=1-E(1,i+1)- E(2,i+1);
end
X=data(320:486,1)';%80-482 3.1-7.1
rate=data(320:486,1)'%
rate(isnan(rate))=0;
n=length(rate);
rt=0:1:n-1;
figure(6)
h11 = plot(rt,X/300000000,'*'); %画点
hold on;
a=E(1,:);
b=E(2,:);
c=E(3,:);
t=1:151;
h1=plot(t,a,'r')%感染人数
%plot(t,b,'b')%健康人数
plot(t,c,'r')%移除者
%legend('i(t)','s(t)','r(t)')
hold on
%% 接种疫苗后
E=zeros(3,150);
l=0.0018; %日接触率 0.0018
m=0.01; %日治愈率 0.01
E(1,1)=data(320,1)/300000000;
E(3,1)=53714928/300000000%43714928 33714928 28664448
E(2,1)=1-E(1,1)-E(3,1);
for i=1:150
E(1,i+1)=E(1,i)+l*E(2,i)-m*E(1,i);
E(2,i+1)=E(2,i)-l*E(1,i)*E(2,i);
E(3,i+1)=1-E(1,i+1)- E(2,i+1);
end
X=data(320:486,1)';%80-482 3.1-7.1
rate=data(320:486,1)'%
rate(isnan(rate))=0;
n=length(rate);
rt=0:1:n-1;
figure(6)
%h11 = plot(rt,X/300000000,'*'); %画点
hold on;
a=E(1,:);
b=E(2,:);
c=E(3,:);
t=1:151;
h2=plot(t,a,'b')%感染人数
%plot(t,b,'b')%健康人数
plot(t,c,'b')%移除者
%legend('i(t)','s(t)','r(t)')
hold on
xlabel('日期'); %设置横坐标名
ylabel('r(t)与i（t）'); %设置纵坐标名
legend([h1 h2 h11],'接种疫苗前','接种疫苗后','感染人数原始数据','Location','NorthWest');
|
