GNU Octave（Matlab）COVID-19 SIR-X模型

Question

我正在尝试执行与本文所述相同的操作：https://science.sciencemag.org/content/sci/early/2020/04/07/science.abb4557.full.pdf

我使用GNU Octave。

这是ODE系统的功能文件：

function f = COVID19ODE(t,y0)

  alpha = 3.07*0.38;
  beta = 0.38;
  kappa = 0.5;
  kappa0 = 0.5;
  ydot =@(t,y,a) ([-alpha*y(1)*y(2) - kappa*y(1);
                  alpha*y(1)*y(2) - beta*y(2) - kappa0*y(2) - kappa*y(2);
                  (kappa0 + kappa)*y(2);
                  kappa0*y(1) + beta*y(2)]);
  odeopt = odeset ("InitialStep", 1e-2, "MaxStep", 1);           
  [t,y] = ode45(ydot, t, y0, odeopt);
  f = y;
endfunction

这是使数据适合ODE模型的脚本：

pkg load optim

Cases = [2,3,20,79,150,227,320,445,650,888,1128,1694,2036,2502,3089,3858,4636,5883,7375,9172,10149,12462,15113,17660,21157,24747,27980,31506,35713,41035,47021,53578,59138,63927,69176,74386,80539,86498,92472,97689,101739,105792,110574,115242,119827,124632,128948,132547,135586,139422,143626,147577,152271];
Days = (1:1:length(Cases));

% Fit
xdata = Days;
ydata = Cases';
F = @(a,xdata) COVID19ODE(xdata,a)(:,3);
a0=[1 1 1 1];
[a,resnorm,~,exitflag,output] = lsqcurvefit(F,a0,xdata,ydata);

% Plot
plot(Days,Cases,'-+' , Days,F(a,Days)) 
grid
legend("Cases" , "Fit" , "location","northwest")
title('Covid 19 pandemic')
xlabel('Days')
ylabel('Cases')

计算时间很长，最后出现错误消息：

warning:  Solving was not successful.  The iterative integration loop exited at time t = 1.000000 before the endpoint at tend= 53.000000 was reached.  This may happen if the stepsize becomes too small.  Try to reduce the value of 'InitialStep' and/or'MaxStep' with the command 'odeset'.
warning: called from
    integrate_adaptive at line 312 column 7
    ode45 at line 232 column 12
    COVID19ODE at line 12 column 8
    COVID19Model>@<anonymous> at line 9 column 16
    nonlin_curvefit>@<anonymous> at line 84 column 14
    __nonlin_residmin__>@<anonymous> at line 316 column 41
    __lm_svd__ at line 469 column 11
    __nonlin_residmin__ at line 452 column 21
    nonlin_curvefit at line 83 column 18
    lsqcurvefit at line 268 column 19
    COVID19Model at line 11 column 30

我试图减小'InitialStep'和/或'MaxStep'的值，但这没有帮助。拟合函数对于数据似乎不是一个很好的近似值：Plot

Answer 1

总结本文，在密度加起来等于S+I+X+R=1的情况下，模型具有模型的“最佳”系数，该系数等于初始值的3个自由参数。在IC中假设X=R=0时，保留一个参数，因此可以将初始条件设置为S0 = 0.9+0.1/(1+x^2)和I=1-S0。此外，假定案例编号与X分量成比例，这将为比例因子提供第二个参数。

通过kappa = 0.2/(1+5*exp(u))使kappas变量在0到0.2的范围内。即使这样，kappa0的使用也会被拒绝，并使用参数

prop =  49040365.30572
y0 =

   1.0000e+00   1.1947e-06   1.0000e-16   1.0000e-16

kappa =  0.042094
kappa0 = 0

我得到对数图

[在参数化的许多其他变体中，求解器找不到任何明智的解，通常，如果根本找不到解，则S分量下降得太快。

此图像的代码是

  % data and other initializations left out
  C = @(Y) Y(:,3);
  L = @(Y) log(1+Y);
  F = @(a,xdata) L(C(COVID19ODE(xdata,a)));

  a0=[0.1 log(1e9) 0 0 ];
  [a,resnorm,~,exitflag,output] = lsqcurvefit(F,a0,xdata,L(ydata));
  disp("a="); disp(a);
  % Plot
  time=Days(1):0.1:Days(end);
  semilogy(Days,Cases,'+r' , time,COVID19ODE(time,a)); 
  grid
  legend("Cases" , "S","I","X","R" ); % , "location","northwest")
  title('Covid 19 pandemic')
  xlabel('Days')
  ylabel('Cases')

  function f = COVID19ODE(t,para)
    ic = para(1);
    prop = exp(para(2))
    S0 =  0.9+0.1/(1+ic^2); 
    y0 = [ S0-2e-16, 1-S0, 1e-16, 1e-16 ]
    beta = 0.38;
    alpha = 3.07*beta;
    kappa = 0.2/(1+5*exp(para(3)))
    kappa0 = 0.2/(1+5*exp(para(4)))
    ydot =@(t,y,a) ([-alpha*y(1)*y(2) - kappa*y(1);
                    alpha*y(1)*y(2) - beta*y(2) - kappa0*y(2) - kappa*y(2);
                    (kappa0 + kappa)*y(2);
                    kappa0*y(1) + beta*y(2)]);
    odeopt = odeset ("InitialStep", 1e-2, "MaxStep", 0.1, "AbsTol", 1e6);           
    [t,y] = ode45(ydot, t, y0, odeopt);
    f = prop*y;
  end%function