其实这个代码四天前就写好了打算发布了,怕数据误差太大打脸,观察了几天数据,发现非常准确,发布。
数学模型是我自己研发的(应该没有人比我先吧,闭门造车不知道外面的情况)
SEIRD
模型简介:
S
表示易感者,即还没感染的人E
表示传染者,即携带有病毒可以传染但是自己没有症状I
表示感染者,即自己感染能感染别人的人R
表示康复者,即已经恢复不会再感染了的人D
表示死亡者,即已经死亡不会传染别人的人beta2
表示传染因子(可以理解为概率)r
表示每个人接触的人数beta1
表示传染者变成感染者的概率gamma
表示感染者康复的概率alpha
表示感染者死亡的概率
然后拟合就好啦
拟合代码:
const days=19;
var
//sti:array[1..days]of longint=(830,1287,1985,2761,4535,5997,7736,9720,11821,14411);
std:array[1..days]of longint=(1,9,17,25,41,56,80,106,132,170,213,259,304,361,425,491,564,637,723);
//str:array[1..days]of longint=(0,0,2,34,38,49,51,60,103,126,171,243,328,475,632,892);
s,e,i,r,d:array[1..days]of extended;
x1,x2,x3,x4,x5,x6,xx1,xx2,xx3,xx4,xx5,xx6,temp,nowmin,testmin:extended;
//j:longint;
tot:longint=0;
function _abs(a,b:extended):extended;begin exit(sqr(a-b){/sqr(b)});end;
function test(alpha,beta1,beta2r,gamma,S0,E0:extended;_out:boolean):extended;
var
now,j:longint;
ans:extended;
begin
s[1]:=S0;e[1]:=E0;i[1]:=830{sti[1]};r[1]:=0{str[1]};d[1]:=std[1];
for now:=2 to days do
begin
s[now]:=s[now-1]-(s[now-1]*(i[now-1]+e[now-1]*beta2r))/(s[now-1]+i[now-1]+e[now-1]);
e[now]:=e[now-1]+(s[now-1]*(i[now-1]+e[now-1]*beta2r))/(s[now-1]+i[now-1]+e[now-1])-e[now-1]*beta1;
i[now]:=i[now-1]+e[now-1]*beta1-gamma*i[now-1]-alpha*i[now-1];
r[now]:=r[now-1]+gamma*i[now-1];
d[now]:=d[now-1]+alpha*i[now-1];
end;
ans:=0;
for j:=1 to days do
begin
//ans:=ans+_abs(i[j],sti[j]);
ans:=ans+_abs(d[j],std[j]);
//ans:=ans+_abs(r[j],str[j]);
end;
if _out then
for j:=1 to days do
writeln({i[j]:20:10,' ',}d[j]:20:10{,' ',r[j]:20:10});
exit(ans);
end;
function _ok(a:extended):boolean;
begin
//writeln('_ok:',a:0:10);
if a>1 then
exit(false);
//exit(random>exp(a-1));
exit(false);
end;
begin
randomize;
x1:=0.00437225297333588;
x2:=0.00824671967809310;
x3:=0.00007918320441526;
x4:=0.00017262652917990;
x5:=293927599.43554592000000000;
x6:=105804.63287489372000000;
xx1:=x1;xx2:=x2;xx3:=x3;xx4:=x4;xx5:=x5;xx6:=x6;
nowmin:=test(x1,x2,x3,x4,x5,x6,true);
//for j:=1 to 100000000 do
//j:=0;
temp:=0.5;
while true do
begin
//j:=j+1;
temp:=temp*0.999998;
//writeln('temp: ',temp:0:10);
{writeln('debug: options: ',xx1:0:20);
writeln('debug: options: ',xx2:0:20);
writeln('debug: options: ',xx3:0:20);
writeln('debug: options: ',xx4:0:20);
writeln('debug: options: ',xx5:0:20);
writeln('debug: options: ',xx6:0:20);
writeln('-------------------------b');}
case random(15) of
0:xx1:=x1*((random*2-1)*temp+1);
1,6,7,8:xx2:=x2*((random*2-1)*temp+1);
2,9,10,11:xx3:=x3*((random*2-1)*temp+1);
3,12,13,14:xx4:=x4*((random*2-1)*temp+1);
4:xx5:=x5*((random*2-1)*temp+1);
5:xx6:=x6*((random*2-1)*temp+1);
end;
{writeln('debug: options: ',xx1:0:20);
writeln('debug: options: ',xx2:0:20);
writeln('debug: options: ',xx3:0:20);
writeln('debug: options: ',xx4:0:20);
writeln('debug: options: ',xx5:0:20);
writeln('debug: options: ',xx6:0:20);
writeln('-------------------------e');}
if xx5>1000000000 then
xx5:=700000000;
testmin:=test(xx1,xx2,xx3,xx4,xx5,xx6,false);
if testmin<nowmin then
begin
//writeln('maj entered');
//readln;
x1:=xx1;x2:=xx2;x3:=xx3;x4:=xx4;x5:=xx5;x6:=xx6;
tot:=tot+1;nowmin:=testmin;
if tot mod 100000=0 then
begin
nowmin:=test(x1,x2,x3,x4,x5,x6,true);
writeln('nowmin = ',nowmin:0:20);
writeln('options: ',x1:0:20);
writeln('options: ',x2:0:20);
writeln('options: ',x3:0:20);
writeln('options: ',x4:0:20);
writeln('options: ',x5:0:20);
writeln('options: ',x6:0:20);
end;
end
else
if _ok((testmin-nowmin)/nowmin) then
begin
writeln('sub entered');
//readln;
x1:=xx1;x2:=xx2;x3:=xx3;x4:=xx4;x5:=xx5;x6:=xx6;
nowmin:=test(x1,x2,x3,x4,x5,x6,true);
writeln('nowmin = ',nowmin:0:20);
writeln('options: ',x1:0:20);
writeln('options: ',x2:0:20);
writeln('options: ',x3:0:20);
writeln('options: ',x4:0:20);
writeln('options: ',x5:0:20);
writeln('options: ',x6:0:20);
end;
end;
end.
然后把输出的几个option在拟合地差不多了(nowmin是误差)时候可以拷到下一份代码里进行未来预估:
const days=26;delta=10;
var
s,e,i,r,d:array[0..days+delta]of extended;
x1,x2,x3,x4,x5,x6:extended;
//j:longint;
tot:longint=0;
function test(alpha,beta1,beta2r,gamma,S0,E0:extended;_out:boolean):extended;
var
now,j:longint;
begin
s[1]:=S0;e[1]:=E0;i[1]:=830;r[1]:=0;d[1]:=1;
for now:=2 to days do
begin
s[now]:=s[now-1]-(s[now-1]*(i[now-1]+e[now-1]*beta2r))/(s[now-1]+i[now-1]+e[now-1]);
e[now]:=e[now-1]+(s[now-1]*(i[now-1]+e[now-1]*beta2r))/(s[now-1]+i[now-1]+e[now-1])-e[now-1]*beta1;
i[now]:=i[now-1]+e[now-1]*beta1-gamma*i[now-1]-alpha*i[now-1];
r[now]:=r[now-1]+gamma*i[now-1];
d[now]:=d[now-1]+alpha*i[now-1];
end;
if _out then
for j:=1 to days do
writeln(d[j]:17:10{,' ',r[j]:17:10,' ',i[j]:17:10}{,' ',s[j]:12:10,' ',e[j]:12:10});
end;
begin
randomize;
x1:=0.00412843138313884 ;
x2:=0.00725008506802688 ;
x3:=0.00010734449832662 ;
x4:=0.00006106562274054 ;
x5:=302215922.84959049000000000;
x6:=129753.99907592588000000 ;
test(x1,x2,x3,x4,x5,x6,true);
end.
以及是否把感染人数等一些可能统计不完全的数据一起拟合进来(这样的话死亡人数预估可能不太准)
数据仅代表模型预测,不代表作者主观观点,不可用于商业用途