SI，SIS，SIR模型的正确实现（python）

Question

我已经创建了上述模型的一些非常基本的实现。但是，尽管图表看起来看起来很正确，但这些数字并不等于常数。这是因为每个隔室中易感染/感染/恢复的人的总和应该总计为N（这是人的总数），但是由于某些原因，它加起来有一些奇怪的十进制数，我真的不知道如何解决它，现在看了3天后。

SI型：

import matplotlib.pyplot as plt 

N = 1000000 
S = N - 1 
I = 1 
beta = 0.6 

sus = [] # infected compartment 
inf = [] # susceptible compartment 
prob = [] # probability of infection at time t 

def infection(S, I, N): 
    t = 0 
    while (t < 100): 
     S = S - beta * ((S * I/N)) 
     I = I + beta * ((S * I)/N) 
     p = beta * (I/N) 

     sus.append(S) 
     inf.append(I) 
     prob.append(p) 
     t = t + 1 

infection(S, I, N) 
figure = plt.figure() 
figure.canvas.set_window_title('SI model') 

figure.add_subplot(211) 
inf_line, =plt.plot(inf, label='I(t)') 

sus_line, = plt.plot(sus, label='S(t)') 
plt.legend(handles=[inf_line, sus_line]) 

plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0)) # use scientific notation 

ax = figure.add_subplot(212) 
prob_line = plt.plot(prob, label='p(t)') 
plt.legend(handles=prob_line) 

type(ax) # matplotlib.axes._subplots.AxesSubplot 

# manipulate 
vals = ax.get_yticks() 
ax.set_yticklabels(['{:3.2f}%'.format(x*100) for x in vals]) 

plt.xlabel('T') 
plt.ylabel('p') 

plt.show() 

SIS型号：

import matplotlib.pylab as plt 

N = 1000000 
S = N - 1 
I = 1 
beta = 0.3 
gamma = 0.1 

sus = \[\] 
inf = \[\] 

def infection(S, I, N): 
    for t in range (0, 1000): 
     S = S - (beta*S*I/N) + gamma * I 
     I = I + (beta*S*I/N) - gamma * I 

     sus.append(S) 
     inf.append(I) 


infection(S, I, N) 

figure = plt.figure() 
figure.canvas.set_window_title('SIS model') 

inf_line, =plt.plot(inf, label='I(t)') 

sus_line, = plt.plot(sus, label='S(t)') 
plt.legend(handles=\[inf_line, sus_line\]) 

plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0)) 

plt.xlabel('T') 
plt.ylabel('N') 

plt.show() 

SIR型号：

import matplotlib.pylab as plt 

N = 1000000 
S = N - 1 
I = 1 
R = 0 
beta = 0.5 
mu = 0.1 

sus = [] 
inf = [] 
rec = [] 

def infection(S, I, R, N): 
    for t in range (1, 100): 
     S = S -(beta * S * I)/N 
     I = I + ((beta * S * I)/N) - R 
     R = mu * I 

     sus.append(S) 
     inf.append(I) 
     rec.append(R) 

infection(S, I, R, N) 

figure = plt.figure() 
figure.canvas.set_window_title('SIR model') 

inf_line, =plt.plot(inf, label='I(t)') 

sus_line, = plt.plot(sus, label='S(t)') 

rec_line, = plt.plot(rec, label='R(t)') 
plt.legend(handles=[inf_line, sus_line, rec_line]) 

plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0)) 

plt.xlabel('T') 
plt.ylabel('N') 


plt.show() 

Answer 1

我只看了SI模型。

您的两个关键变量是S和I。 （你可能会颠倒这两个变量的含义，尽管这并不影响我在这里写的东西。）你将它们初始化，使它们的总和为N，这是常数1000000。

您在线路

S = S - beta * ((S * I/N)) 
I = I + beta * ((S * I)/N) 

你显然打算增加I和S减去相同的值更新你的两个关键变量，所以S和I总和不变。但是，您实际上首先更改S，然后使用该新值更改I，因此添加和减去的值实际上并不相同，并且变量的总和未保持不变。

您可以通过使用Python更新多行变量的能力来解决此问题。与

S, I = S - beta * ((S * I/N)), I + beta * ((S * I)/N) 

这两种计算新值的更新变量替换前这两条线，所以相同的数值实际上添加，并从两个变量中减去。 （还有其他方法可以获得相同的效果，例如更新值的临时变量，或者一个临时变量来存储要添加和减去的数量，但由于您使用Python，因此您可以使用其功能。）

当我现在运行的程序，我得到这些图表：

我认为这是你想要的。

Answer 2

因此，上述解决方案也适用于SIS模型。

至于SIR模型我不得不使用odeint解微分方程，这里是一个简单的解决方案SIR模型：

import matplotlib.pylab as plt 
from scipy.integrate import odeint 
import numpy as np 

N = 1000 
S = N - 1 
I = 1 
R = 0 
beta = 0.6 # infection rate 
gamma = 0.2 # recovery rate 

# differential equatinons 
def diff(sir, t): 
    # sir[0] - S, sir[1] - I, sir[2] - R 
    dsdt = - (beta * sir[0] * sir[1])/N 
    didt = (beta * sir[0] * sir[1])/N - gamma * sir[1] 
    drdt = gamma * sir[1] 
    print (dsdt + didt + drdt) 
    dsirdt = [dsdt, didt, drdt] 
    return dsirdt 


# initial conditions 
sir0 = (S, I, R) 

# time points 
t = np.linspace(0, 100) 

# solve ODE 
# the parameters are, the equations, initial conditions, 
# and time steps (between 0 and 100) 
sir = odeint(diff, sir0, t) 

plt.plot(t, sir[:, 0], label='S(t)') 
plt.plot(t, sir[:, 1], label='I(t)') 
plt.plot(t, sir[:, 2], label='R(t)') 

plt.legend() 

plt.xlabel('T') 
plt.ylabel('N') 

# use scientific notation 
plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0)) 

plt.show()