正确实施 SI、SIS、SIR 模型 (python)
Correct implementation of SI, SIS, SIR models (python)
我已经创建了上述模型的一些非常基本的实现。然而，尽管图表看起来是正确的，但数字加起来并不是一个常数。那是因为每个车厢里 susceptible/infected/recovered 人的总和应该加起来为 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()
我只看 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 你也可以使用它的功能。 )
当我现在 运行 程序时，我得到这些图表：
我想这就是你想要的。
所以上面的解决方案也适用于 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()
我已经创建了上述模型的一些非常基本的实现。然而，尽管图表看起来是正确的，但数字加起来并不是一个常数。那是因为每个车厢里 susceptible/infected/recovered 人的总和应该加起来为 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()
我只看 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 你也可以使用它的功能。 )
当我现在 运行 程序时，我得到这些图表：
我想这就是你想要的。
所以上面的解决方案也适用于 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()