受杰森的《Almost Looks Like Work》啟發(fā)，我來展示一些病毒傳播模型。需要注意的是這個模型并不反映現(xiàn)實情況，因此不要誤以為是西非可怕的傳染病。相反，它更應(yīng)該被看做是某種虛構(gòu)的僵尸爆發(fā)現(xiàn)象。那么，讓我們進(jìn)入主題。

2015331152332590.jpg (565×105)

這就是SIR模型，其中字母S、I和R反映的是在僵尸疫情中，個體可能處于的不同狀態(tài)。

S 代表易感群體，即健康個體中潛在的可能轉(zhuǎn)變的數(shù)量。
I 代表染病群體，即僵尸數(shù)量。
R 代表移除量，即因死亡而退出游戲的僵尸數(shù)量，或者感染后又轉(zhuǎn)回人類的數(shù)量。但對與僵尸不存在治愈者，所以我們就不要自我愚弄了（如果要把SIR模型應(yīng)用到流感傳染中，還是有治愈者的）。
至于β(beta)和γ(gamma):
β(beta)表示疾病的傳染性程度，只要被咬就會感染。
γ(gamma)表示從僵尸走向死亡的速率，取決于僵尸獵人的平均工作速率，當(dāng)然，這不是一個完美的模型，請對我保持耐心。
S′=?βIS告訴我們健康者變成僵尸的速率，S′是對時間的導(dǎo)數(shù)。
I′=βIS?γI告訴我們感染者是如何增加的，以及行尸進(jìn)入移除態(tài)速率（雙關(guān)語）。
R′=γI只是加上(gamma I)，這一項在前面的等式中是負(fù)的。

上面的模型沒有考慮S/I/R的空間分布，下面來修正一下！

一種方法是把瑞典和北歐國家分割成網(wǎng)格，每個單元可以感染鄰近單元，描述如下：

其中對于單元，和是它周圍的四個單元。（不要因為對角單元而腦疲勞，我們需要我們的大腦不被吃掉）。

初始化一些東東。

import numpy as np
import math
import matplotlib.pyplot as plt  
%matplotlib inline
from matplotlib import rcParams
import matplotlib.image as mpimg
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16
rcParams['figure.figsize'] = 12, 8
from PIL import Image

適當(dāng)?shù)腷eta和gamma值就能夠摧毀大半江山

beta = 0.010
gamma = 1

還記得導(dǎo)數(shù)的定義么？當(dāng)導(dǎo)數(shù)已知，假設(shè)Δt很小的情況下，經(jīng)過重新整理，它可以用來近似預(yù)測函數(shù)的下一個取值，我們已經(jīng)聲明過u′(t)。

2015331152448897.jpg (631×171)

初始化一些東東。

import numpy as np
import math
import matplotlib.pyplot as plt  
%matplotlib inline
from matplotlib import rcParams
import matplotlib.image as mpimg
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16
rcParams['figure.figsize'] = 12, 8
from PIL import Image

適當(dāng)?shù)腷eta和gamma值就能夠摧毀大半江山

beta = 0.010
gamma = 1

還記得導(dǎo)數(shù)的定義么？當(dāng)導(dǎo)數(shù)已知，假設(shè)Δt很小的情況下，經(jīng)過重新整理，它可以用來近似預(yù)測函數(shù)的下一個取值，我們已經(jīng)聲明過u′(t)。

2015331152522306.jpg (549×363)

這種方法叫做歐拉法，代碼如下：

def euler_step(u, f, dt):
  return u + dt * f(u)

我們需要函數(shù)f(u)。友好的numpy提供了簡潔的數(shù)組操作。我可能會在另一篇文章中回顧它，因為它們太強(qiáng)大了，需要更多的解釋，但現(xiàn)在這樣就能達(dá)到效果：

 
def f(u):
  S = u[0]
  I = u[1]
  R = u[2]
 
  new = np.array([-beta*(S[1:-1, 1:-1]*I[1:-1, 1:-1] +
              S[0:-2, 1:-1]*I[0:-2, 1:-1] +
              S[2:, 1:-1]*I[2:, 1:-1] +
              S[1:-1, 0:-2]*I[1:-1, 0:-2] +
              S[1:-1, 2:]*I[1:-1, 2:]),
           beta*(S[1:-1, 1:-1]*I[1:-1, 1:-1] +
              S[0:-2, 1:-1]*I[0:-2, 1:-1] +
              S[2:, 1:-1]*I[2:, 1:-1] +
              S[1:-1, 0:-2]*I[1:-1, 0:-2] +
              S[1:-1, 2:]*I[1:-1, 2:]) - gamma*I[1:-1, 1:-1],
           gamma*I[1:-1, 1:-1]
          ])
 
  padding = np.zeros_like(u)
  padding[:,1:-1,1:-1] = new
  padding[0][padding[0] < 0] = 0
  padding[0][padding[0] > 255] = 255
  padding[1][padding[1] < 0] = 0
  padding[1][padding[1] > 255] = 255
  padding[2][padding[2] < 0] = 0
  padding[2][padding[2] > 255] = 255
 
  return padding

導(dǎo)入北歐國家的人口密度圖并進(jìn)行下采樣，以便較快地得到結(jié)果

from PIL import Image
img = Image.open('popdens2.png')
img = img.resize((img.size[0]/2,img.size[1]/2))
img = 255 - np.asarray(img)
imgplot = plt.imshow(img)
imgplot.set_interpolation('nearest')

2015331152603003.jpg (373×489)

北歐國家的人口密度圖（未包含丹麥）

S矩陣，也就是易感個體，應(yīng)該近似于人口密度。感染者初始值是0，我們把斯德哥爾摩作為第一感染源。

S_0 = img[:,:,1]
I_0 = np.zeros_like(S_0)
I_0[309,170] = 1 # patient zero

因為還沒人死亡，所以把矩陣也置為0.

R_0 = np.zeros_like(S_0)

接著初始化模擬時長等。

T = 900             # final time
dt = 1             # time increment
N = int(T/dt) + 1        # number of time-steps
t = np.linspace(0.0, T, N)   # time discretization
 
# initialize the array containing the solution for each time-step
u = np.empty((N, 3, S_0.shape[0], S_0.shape[1]))
u[0][0] = S_0
u[0][1] = I_0
u[0][2] = R_0

我們需要自定義一個顏色表，這樣才能將感染矩陣顯示在地圖上。

import matplotlib.cm as cm
theCM = cm.get_cmap("Reds")
theCM._init()
alphas = np.abs(np.linspace(0, 1, theCM.N))
theCM._lut[:-3,-1] = alphas

下面坐下來欣賞吧…

 
for n in range(N-1):
  u[n+1] = euler_step(u[n], f, dt)

讓我們再做一下圖像渲染，把它做成gif，每個人都喜歡gifs！

from images2gif import writeGif
 
keyFrames = []
frames = 60.0
 
for i in range(0, N-1, int(N/frames)):
  imgplot = plt.imshow(img, vmin=0, vmax=255)
  imgplot.set_interpolation("nearest")
  imgplot = plt.imshow(u[i][1], vmin=0, cmap=theCM)
  imgplot.set_interpolation("nearest")
  filename = "outbreak" + str(i) + ".png"
  plt.savefig(filename)
  keyFrames.append(filename)
 
images = [Image.open(fn) for fn in keyFrames]
gifFilename = "outbreak.gif"
writeGif(gifFilename, images, duration=0.3)
plt.clf()

您可能感興趣的文章: