Never too old to learn.

SARI病毒的数据拟合

Posted on By Andy Zhu

其实这个代码四天前就写好了打算发布了,怕数据误差太大打脸,观察了几天数据,发现非常准确,发布。

数学模型是我自己研发的(应该没有人比我先吧,闭门造车不知道外面的情况)

1.png

SEIRD模型简介:

  1. S表示易感者,即还没感染的人
  2. E表示传染者,即携带有病毒可以传染但是自己没有症状
  3. I表示感染者,即自己感染能感染别人的人
  4. R表示康复者,即已经恢复不会再感染了的人
  5. D表示死亡者,即已经死亡不会传染别人的人
  6. beta2表示传染因子(可以理解为概率)
  7. r表示每个人接触的人数
  8. beta1表示传染者变成感染者的概率
  9. gamma表示感染者康复的概率
  10. 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.

以及是否把感染人数等一些可能统计不完全的数据一起拟合进来(这样的话死亡人数预估可能不太准)

数据仅代表模型预测,不代表作者主观观点,不可用于商业用途