目前整個世界同新冠肺炎的鬥爭仍在繼續，從新冠肺炎爆發以來，很多人都提出了各種各樣的模型來預測肺炎疫情的發展，其中比較常見的就是SIR模型。這是因為SIR是目前疾病防控領域最經典也是最常用的一個模型，而今天筆者就用圖論來講述一下SIR模型。
首先對SIR模型做一個簡單的介紹，SIR全稱就是Susceptible—Infected—Recovered，翻譯過來就是易感—感染—康復，即易感人群Susceptible有α概率被某種疾病感染，成為感染人群Infected，而感染人群Infected又有β概率康復，成為康復人群Recovered。這個過程可以是一次或者是多次的，而且還衍生出不少類似的模型，比如SIS、SIRS和SI等。整個SIR過程可以用下圖來描述。
圖1. SIR模型示意圖
而筆者這次用的是圖論，也就是我們平常所說的網絡分析來說明SIR模型，大概過程如下圖所示。
圖2. SIR網絡模型示意圖
這個過程其實很簡單，就是由少數人感染多數人，然後最終所有感染人群都康復的一個過程，這個模型簡化了我們很多實際情況所遇到的部分問題，比如不考慮感染人群的死亡問題，不考慮感染人群的再感染等等。這次用到的圖論方面的庫是networkx，這是Python中最常用的圖論分析的庫。下面還是直接用代碼來具體說明一下。
首先導入各種包。
<code>import matplotlib.pyplot as plt
import numpy.random as rdm
import networkx as nx/<code>
接下來定義兩個變量n和g。n是總人數，而nx.erdos_renyi_graph則是一種隨機網絡圖，Erdos和Renyi是兩位匈牙利數學家，他們的主要貢獻就是建立了著名的E-R隨機圖理論（Random graph theory），這被公認為在數學上開創了複雜網絡拓撲結構的系統性分析。而這個E-R圖也是本文應用的基礎。而nx.erdos_renyi_graph(n, 0.01)的意思就是以概率0.01來連接100個節點。這裡100和0.01這兩個數字筆者用的比較隨意，大家可以根據自己的研究來任意設定相關數值，這裡主要還是用於簡單說明。然後是三個string類型的變量susceptible、infected和recovered，分別代表“易感”、“感染”和“康復”人群。
<code>n = 100
g = nx.erdos_renyi_graph(n, 0.01)
susceptible = 'S'
infected = 'I'
recovered = 'R'/<code>
接下來我們要定義幾個函數。
<code>def onset(g): #初始設置
for i in g.node.keys():
g.node[i]['state'] = susceptible
def infect_prop(g, proportion): #設置感染比例
for i in g.node.keys():
if(rdm.random() <= proportion):
g.node[i]['state'] = infected
def build_model(pInfect, pRecover): #模型構建
def model(g, i):
if g.node[i]['state'] == infected:
for m in g.neighbors(i):
if g.node[m]['state'] == susceptible:
if rdm.random() <= pInfect:
g.node[m]['state'] = infected
if rdm.random() <= pRecover:
g.node[i]['state'] = recovered
return model
def model_run(g, model): #單次模型運行
for i in g.node.keys():
model(g, i)
def model_iter(g, model, iter): #多次模型循環
for i in range(iter):
model_run(g, model)/<code>
第一個函數onset，也就是剛開始的狀態，在這裡我們把每個節點的狀態都設置為“S”，也就是易感的意思，這裡我們假設所有人都屬於易感人群中的。然後就開始有部分人群出現感染，而這個感染函數就是infect_prop，prop意思就是proportion，這個函數的意思就是感染比例，這裡我們用到了numpy.random.random方法，也就是返回一個均勻分佈，數值大小在0-1之間，不包含1，當rdm.random()<= proportion時，我們就讓這樣的人群變為感染人群，這樣讓感染人群更加均勻一些。第三個函數是build_model，是一個嵌套函數，有兩個參數，pInfect和pRecover，分別代表感染概率和康復概率，而這裡的g.neighbors(i)的意思是節點i的每個相鄰的節點，而g是我們生成的一個圖（graph），是一個class，後面會有說明。在這些感染人群中，當他們的相鄰節點的狀態是“susceptible”時，我們讓這些節點（人群）以概率pInfect來進行感染；而在這些感染人群中，以pRecover的概率進行康復。這就是一個簡單的感染—康復過程。而後面的兩個函數model_run和model_iter則是將這個模型運行一次和多次，分別用來模擬一個循環和多個循環。
接下來就是畫圖。
<code>fig, ax= plt.subplots(figsize=(12, 10))
ax.set_xticks([])
ax.set_yticks([])
pos = nx.spring_layout(g, k=0.2)
nx.draw_networkx_edges(g, pos, alpha=0.5, width = 1)
nx.draw_networkx_nodes(g, pos, node_size=80)
plt.show()/<code>
首先設置圖片的大小，並去掉座標軸，然後設置網絡圖的位置，nx.spring_layout就是設置網絡圖的位置的方法，k是節點間的最佳距離，這個可以隨意設置，值越大節點越分散，接下來繪製節點和連線，nx.draw_networkx_nodes(g, pos, node_size=80)用來繪製節點，pos就是剛才設置的位置參數，再設置一下節點的大小，而nx.draw_networkx_edges(g, pos, alpha=0.5, width = 1)用來繪製連線，同樣要傳入位置參數，再設置透明度和線寬。最後生成的圖如下。
圖3. 模型生成的網絡圖
networkx使用的繪圖算法是隨機的，同時我們使用的參數也是隨機的，所以這個圖每次生成的結果都不同，但大體相似，我們可以看到這裡面已經有部分節點相連，疾病也就是通過他們開始傳播。
最後就是計算感染率。
<code>onset(g)
infect_prop(g, 0.05)
model = build_model(0.2, 0.8)
model_iter(g, model, 10)
infected = [ v for (v, attr) in g.nodes(data = True) if attr['state'] == recovered ]
infection_rate = len(infected)/n
print(infection_rate)/<code>
這當中onset(g)中的g就是前面我們用g = nx.erdos_renyi_graph(n, 0.01)生成的圖，這是一個類的實例，infect_prop(g, 0.05)中的0.05就是設置人群初始感染率為0.05，model = build_model(0.2, 0.8)中的感染幾率和康復幾率分別設置為0.2和0.8，至於為何這麼設置，主要是根據常用的“二八理論”，這個數字可以根據使用者的模擬情況隨意設置，而model_iter(g, model, 10)中我們讓這個模擬過程重複10次，最後計算出感染人數，並得出最終穩定的感染率infection_rate，因為這些參數不少都是隨機的，所以得出的結果理論上每次都是不同的，筆者得出的結果從0.03到0.12不等，這個結果意義不大，主要還是理解模擬的方法。