R語言SIR模型(Susceptible Infected Recovered Model)代碼sir模型實(shí)例
原文鏈接:http://tecdat.cn/?p=14593?
原文出處:拓端數(shù)據(jù)部落公眾號(hào)
?
SIR模型定義
SIR模型是一種傳播模型,是信息傳播過程的抽象描述。
SIR模型是傳染病模型中最經(jīng)典的模型,其中S表示易感者,I表示感染者,R表示移除者。
S:Susceptible,易感者
I:Infective,感染者
R:Removal,移除者
SIR模型的應(yīng)用
SIR模型應(yīng)用于信息傳播的研究。
傳播過程大致如下:最初,所有的節(jié)點(diǎn)都處于易感染狀態(tài)。然后,部分節(jié)點(diǎn)接觸到信息后,變成感染狀態(tài),這些感染狀態(tài)的節(jié)點(diǎn)試著去感染其他易感染狀態(tài)的節(jié)點(diǎn),或者進(jìn)入恢復(fù)狀態(tài)。感染一個(gè)節(jié)點(diǎn)即傳遞信息或者對(duì)某事的態(tài)度。恢復(fù)狀態(tài),即免疫,處于恢復(fù)狀態(tài)的節(jié)點(diǎn)不再參與信息的傳播。
SIR的微分方程
a為感染率、b恢復(fù)率
注意:
t為某個(gè)時(shí)刻,例如t=1,S(1)為第一天易感人群的人數(shù)。
無論t為什么時(shí)刻,總?cè)藬?shù)是不變的,即N(t)=S(t)+I(t)+R(t)。
人口總數(shù)總保持一個(gè)常數(shù),即N(t)=K,不考慮人口的出生、死亡、遷移等因素。
這里介紹一個(gè)使用R模擬網(wǎng)絡(luò)擴(kuò)散的例子。
?
第一步,生成網(wǎng)絡(luò)。
規(guī)則網(wǎng)
g =graph.tree(size, children =2); plot(g)
g =graph.star(size); plot(g)
g =graph.full(size); plot(g)
g =graph.ring(size); plot(g)
第二步,隨機(jī)選取一個(gè)或n個(gè)隨機(jī)種子。
?
# initiate the diffusers
seeds_num =1 diffusers =sample(V(g),seeds_num) ;
diffusers
## + 1/50 vertex:
## [1] 43
infected =list()
infected[[1]]=diffusers#
?
第三步,傳染能力
在這個(gè)簡(jiǎn)單的例子中,每個(gè)節(jié)點(diǎn)的傳染能力是0.5,即與其相連的節(jié)點(diǎn)以0.5的概率被其感染,每個(gè)節(jié)點(diǎn)的回復(fù)能力是0.5,即其以0.5的概率被其回復(fù)。在R中的實(shí)現(xiàn)是通過拋硬幣的方式來實(shí)現(xiàn)的。
## [1] 0
顯然,這很容易擴(kuò)展到更一般的情況,比如節(jié)點(diǎn)的平均感染能力是0.128,那么可以這么寫: 節(jié)點(diǎn)的平均回復(fù)能力是0.1,那么可以這么寫
p =0.128
coins =c(rep(1, p*1000), rep(0,(1-p)*1000))
sample(coins, 1, replace=TRUE, prob=rep(1/n, n))
## [1] 0
n =length(coins2)
sample(coins2, 1, replace=TRUE, prob=rep(1/n, n))
## [1] 0
當(dāng)然最重要的一步是要能按照“時(shí)間”更新網(wǎng)絡(luò)節(jié)點(diǎn)被感染的信息。
?
keep =unlist(lapply(nearest_neighbors[,2], toss))
new_infected =as.numeric(as.character(nearest_neighbors[,1][keep >=1]))
diffusers =unique(c(as.numeric(diffusers), new_infected))
return(diffusers)}
set.seed(1);
?
開啟擴(kuò)散過程!
先看看S曲線吧:
為了可視化這個(gè)擴(kuò)散的過程,我們用紅色來標(biāo)記被感染者。
# generate a palette#
plot(g, layout =layout.old)
set.seed(1)#
library(animation)# start the plot
m =1
如同在Netlogo里一樣,我們可以把網(wǎng)絡(luò)擴(kuò)散與增長(zhǎng)曲線同時(shí)展示出來:
set.seed(1)
# start the plot
m =1
p_cum=numeric(0)
h_cum=numeric(0)
i_cum=numeric(0)
while( m<50 ) {# start the plot
layout(matrix(c(1, 2, 1, 3), 2,2, byrow =TRUE), widths=c(3,1), heights=c(1, 1))
V(g)$color = "white"
V(g)$color[V(g)%in%infected[[m ]] ] = "red"
V(g)$color[V(g)%in%health[[m ]]] = "green"
if(m<=length(infected))
plot(pp~time, type ="h", ylab ="PDF", xlab ="Time",xlim =c(0,i), ylim =c(0,1), frame.plot =FALSE)
m =m +1
}
?
參考文獻(xiàn)
1.R語言泊松Poisson回歸模型分析案例
2.R語言進(jìn)行數(shù)值模擬:模擬泊松回歸模型
3.r語言泊松回歸分析
4.R語言對(duì)布豐投針(蒲豐投針)實(shí)驗(yàn)進(jìn)行模擬和動(dòng)態(tài)可視化
5.用R語言模擬混合制排隊(duì)隨機(jī)服務(wù)排隊(duì)系統(tǒng)
6.GARCH(1,1),MA以及歷史模擬法的VaR比較
7.R語言做復(fù)雜金融產(chǎn)品的幾何布朗運(yùn)動(dòng)的模擬
8.R語言進(jìn)行數(shù)值模擬:模擬泊松回歸模型
9.R語言對(duì)巨災(zāi)風(fēng)險(xiǎn)下的再保險(xiǎn)合同定價(jià)研究案例:廣義線性模型和帕累托分布Pareto distributions