此文用于整理回顾写论文时看的文献资料和学到的知识，也希望能带来一些参考。

什么是SIR模型？

1927年，Kermack 和 McKendrick 为了研究17世纪肆虐伦敦的黑死病和20世纪席卷孟买的鼠疫而发展出的仓室模型（compartmental model）就是现在被广泛应用的SIR模型。其中S代表Susceptible 易感者，I 代表 Infectious 感染者，R代表Recovered/Removed 康复者/移除者。之间的关系可用下图描述，SIR的模型建立在人口不变即无出生无死亡无移民的情况下，并且人口均匀混合在一起（homogeneous mixing。假设总人口为 N N N,那么 N = S + I + R N=S+I+R N=S+I+R.

这三种状态可用常微分方程表示为
d S d t = − β S I \frac{dS}{dt}= - \beta S I dtdS​=−βSI
d I d t = β S I − γ I \frac{dI}{dt}= \beta S I -\gamma I dtdI​=βSI−γI
d R d t = γ S I \frac{dR}{dt}= \gamma S I dtdR​=γSI

其中， β \beta β是传播率，这也暗示感染者仅会和易感者接触，传染率越大，传播越快。 γ \gamma γ是康复率， 1 γ \frac{1}{\gamma} γ1​ 则是感染者携带病毒的平均时间段。新冠病毒通常为5-14天。

其他重要概念还有传染力 λ \lambda λ 指的是从易感者到感染者的速度。 λ = β I / N \lambda=\beta I/N λ=βI/N。另一个重要概念是基本感染数 R 0 R_{0} R0​，指在没有任何干预和免疫的情况下，一个人会把疾病传染给多少人的平均数。这是在疫情初期很多传染病学专家和学者科普过的概念， R 0 R_{0} R0​是通过模型计算得出， R 0 = β / γ R_{0}=\beta / \gamma R0​=β/γ,所以不同文献中所计算出来的值不一样也是很普遍的，因为不同学者的建模不同。如果 R 0 R_{0} R0​小于1，该传染病不会爆发，如果其大于1，则会在社群中爆发。与 R 0 R_{0} R0​相似但并不相同的另一个重要概念是 R t R_{t} Rt​有效感染数。 R t R_{t} Rt​则是建立在存在一定免疫力的情况下，一个人平均会传染多少人。群体免疫这个概念在初期被英国的公共卫生学家大肆宣扬，但是这要建立在疫苗普及的基础上，通过接种疫苗或感染获得免疫力。缺乏免疫力的脆弱个体被感染的机会免疫力的脆弱个体被感染的机会就会减少。

SIR 模型变体

基于上述假设，自SIR模型建立以来，后来的学者们又提出了一些别的模型，如SI模型、SIS模型、SIRS模型、SEIR模型等。这些模型加速了对传染病的研究，促进了对传染病的认识。SI模型用于难以治愈的疾病，例如艾滋病。SIS模型适用于
适用于可以治愈的传染病，如常见的流感。
流感。SIRS模型适用于有二次感染风险的疾病。
SEIR模型适用于被感染的人被治愈并具有免疫力的情况，通常用于麻疹、腮腺炎、风疹等。

确定性模型和随机性模型

SIR模型又可以分为确定性和随机性模型，SIR随机模型是在SIR决定性模型的基础上构建的。所以这两个模型有相同的假设。因此，随机模型中的表达式
模型中的表达式仍然遵循确定性模型中的表达式。

随机性模型不同在被感染的个体在一个随机的时期内保持感染性，感染期被表示为一个独立的、相同的分布。带有强度参数的指数分布在文献中受到了特别的关注。如果未来的状态取决于现在和过去的状态，那么随机过程就是马尔科夫的或具有马尔科夫（Markov Chain）的特性。除了马尔科夫过程之外
随机过程也可以用非马尔科夫方法来模拟，如Sellke construction。

这两个模型的主要区别在于背后的数学原理；经典模型是由普通微分方程驱动的，而随机模型是由连续时间的马尔可夫链或非马尔科夫过程驱动。尽管SIR确定性模型可以成功地对某些传染病进行建模并获得相对较好的结论。但作为随机模型，它的效果并不理想，曲线过于平滑。随机模型的建立成功填补这一空白，以应对某些具有更多可变因素和限制条件的传染病。

下图为确定性模型中S,I,R的变化。N=1000

R语言应用

library(deSolve)
sir_model<-function(time,state,parameters){
  with(as.list(c(state,parameters)),{
    N<-S+I+R
    lambda=beta*I/N #force of infection
    dS=-lambda*S
    dI=lambda*S-gamma*I
    dR=gamma*I
    return(list(c(dS,dI,dR)))
  })
}

N<-1000
state_values<-c(S=N-1,I=1)
r0<-2
t<-10
parameters<-c(beta=r0/t,gamma=1/t)
#timeframe
times<-seq(1,150,by=1)

output<-as.data.frame(ode(y=state_values,
                          times=t,
                          func=sir_model,
                          parms=parameters))