Модель GNU Octave (Matlab) COVID-19 SIR-X - PullRequest
Новая
       76

Модель GNU Octave (Matlab) COVID-19 SIR-X

0 голосов
/ 12 апреля 2020

Я пытаюсь сделать то же самое, что описано в этой статье: 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, но это не помогло. Функция подбора не очень подходит для данных: График

1 Ответ

0 голосов
/ 16 апреля 2020

Подводя итог статьи, мы имеем «лучшие» коэффициенты модели для случая плотностей, которые складываются в единицу, S+I+X+R=1, что составляет 3 свободных параметра для начального значения. При допущении X=R=0 в I C остается один параметр, так что начальные условия можно параметризировать как S0 = 0.9+0.1/(1+x^2) и I=1-S0. Кроме того, предполагается, что номер дела пропорционален компоненту X, которая даст второй параметр для коэффициента пропорциональности.

Установите переменную Каппа в диапазоне от 0 до 0,2 через kappa = 0.2/(1+5*exp(u)). Даже тогда использование kappa0 отклонено, с параметрами 

prop =  49040365.30572
y0 =

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

kappa =  0.042094
kappa0 = 0

Я получаю логарифмический c участок enter image description here

Во многих других вариантах При параметризации решатель не находит разумного решения, обычно, если решение вообще найдено, компонент 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
...