代码之家 › 专栏 › 技术社区 › Terry BRETT

SI、SIS、SIR模型的正确实现(python)

modeling math python

Terry BRETT · 技术社区 · 7 年前

我已经创建了上述模型的一些非常基本的实现。然而,尽管图表看起来很正确,但这些数字加起来并不是一个常数。也就是说,每个隔间中易感/感染/康复者的总数应该加起来是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()

2 回复 | 直到 7 年前

Rory Daulton 7 年前

我只看SI模型。

你的两个关键变量是 S 和 I . (你可能颠倒了这两个变量的含义,但这并不影响我在这里写的内容。)初始化它们,使其总和为 N 这是常数 1000000 .

您更新了行中的两个关键变量

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

你显然想增加 我 并从中减去相同的值,因此 S 和 我 保持不变。然而,你实际上首先改变了 S 然后使用该新值进行更改 我 ,因此加和减的值实际上并不相同,变量之和也没有保持恒定。

您可以通过使用Python在一行中更新多个变量的能力来解决这个问题。将这两条线替换为

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

这会在更新变量之前计算两个新值,因此实际上从两个变量中添加和减去的值是相同的。(还有其他方法可以获得相同的效果,例如用于更新值的临时变量,或者一个临时变量来存储要加和减的量,但是由于使用Python,您也可以使用它的功能。)

我想这就是你想要的。

Terry BRETT 7 年前

因此,上述解决方案也适用于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()