原文鏈接：http://tecdat.cn/?p=14593?

與普通的擴(kuò)散研究不同，網(wǎng)絡(luò)擴(kuò)散開始考慮網(wǎng)絡(luò)結(jié)構(gòu)對(duì)于擴(kuò)散過程的影響。 這里介紹一個(gè)使用R模擬網(wǎng)絡(luò)擴(kuò)散的例子。

基本的算法非常簡(jiǎn)單： 生成一個(gè)網(wǎng)絡(luò):g(V, E)。 隨機(jī)選擇一個(gè)或幾個(gè)節(jié)點(diǎn)作為種子（seeds）。 每個(gè)感染者以概率p（可視作該節(jié)點(diǎn)的傳染能力,通常表示為ββ）影響與其相連的節(jié)點(diǎn)。 其實(shí)這是一個(gè)最簡(jiǎn)單的SI模型在網(wǎng)絡(luò)中的實(shí)現(xiàn)。S表示可感染（susceptible）, I表示被感染（infected）。易感態(tài)-感染態(tài)-恢復(fù)態(tài)(SIR)模型用以描述水痘和麻疹這類患者能完全康復(fù)并獲得終身免疫力的流行病。對(duì)于SIR流行病傳播模型，任意時(shí)刻節(jié)點(diǎn)只能處于易感態(tài)(S)或感染態(tài)(I)或恢復(fù)態(tài)(R)。易感態(tài)節(jié)點(diǎn)表示未被流行病感染的個(gè)體，且可能被感染；感染態(tài)節(jié)點(diǎn)表示已經(jīng)被流行病感染且具有傳播能力；恢復(fù)態(tài)節(jié)點(diǎn)則表示曾感染流行病且完全康復(fù)。與SIS模型類似，每一時(shí)間步內(nèi)，每個(gè)感染態(tài)節(jié)點(diǎn)以概率λλ嘗試感染它的鄰居易感態(tài)節(jié)點(diǎn)，并以概率γγ變?yōu)榛謴?fù)態(tài)。SIR模型可以表達(dá)為：

S = S（t）是易感個(gè)體的數(shù)量， I = I（t）是被感染的個(gè)體的數(shù)目， R = R（t）是恢復(fù)的個(gè)體的數(shù)目。

第二組因變量代表在三個(gè)類別的總?cè)丝诘谋壤?。所以，如果N是總?cè)丝冢?90萬在我們的例子），我們有

S（T）= S（T）/ N，人口的易感部分， Ⅰ（T）= I（t）的/ N的人口感染分?jǐn)?shù)并 R（T）= R（t）的/ N，人口的康復(fù)部分。

解這個(gè)微分方程，我們可以得到累計(jì)增長(zhǎng)曲線的表達(dá)式。有趣的是，這是一個(gè)logistic增長(zhǎng)，具有明顯的S型曲線（S-shaped curve）特征。該模型在初期跨越臨界點(diǎn)之后增長(zhǎng)較快，后期則變得緩慢。 因而可以用來描述和擬合創(chuàng)新擴(kuò)散過程（diffusion of innovations）。 當(dāng)然，對(duì)疾病傳播而言，SI模型是非常初級(jí)的（naive），主要因?yàn)槭芨腥镜膫€(gè)體以一定的概率恢復(fù)健康，或者繼續(xù)進(jìn)入可以被感染狀態(tài)(S，據(jù)此擴(kuò)展為SIS模型)或者轉(zhuǎn)為免疫狀態(tài)（R,據(jù)此擴(kuò)展為SIR模型）。 免疫表示為R，用γγ代表免疫概率（removal or recovery rate)。對(duì)于信息擴(kuò)散而言，這種考慮暫時(shí)是不需要的。

?

第一步，生成網(wǎng)絡(luò)。

規(guī)則網(wǎng)

1. g =graph.tree(size, children =2); plot(g)

2. 


1. g =graph.star(size); plot(g)

2. 


1. g =graph.full(size); plot(g)

2. 


1. g =graph.ring(size); plot(g)

2. 


1. g =connect.neighborhood(graph.ring(size), 2); plot(g) # 最近鄰耦合網(wǎng)絡(luò)

2. 


1. # 隨機(jī)網(wǎng)絡(luò)g =erdos.renyi.game(size, 0.1)# 小世界網(wǎng)絡(luò)

2. 


3. g = rewire.edges(erdos.renyi.game(size, 0.1), prob = 0.8 )# 無標(biāo)度網(wǎng)絡(luò)

4. 


5. g =barabasi.game(size) ; plot(g)

第二步，隨機(jī)選取一個(gè)或n個(gè)隨機(jī)種子。
?

1. # initiate the diffusers

2. 


3. seeds_num =1 diffusers =sample(V(g),seeds_num) ;

4. 


5. diffusers

6. 


7. ## + 1/50 vertex:

8. 


9. ## [1] 43

10. 


11. infected =list()

12. 


13. 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. ## [1] 0

2. 


顯然，這很容易擴(kuò)展到更一般的情況，比如節(jié)點(diǎn)的平均感染能力是0.128，那么可以這么寫： 節(jié)點(diǎn)的平均回復(fù)能力是0.1，那么可以這么寫

1. p =0.128

2. 


3. coins =c(rep(1, p*1000), rep(0,(1-p)*1000))

4. 


5. sample(coins, 1, replace=TRUE, prob=rep(1/n, n))

6. 


7. ## [1] 0

8. 


9. n =length(coins2)

10. 


11. sample(coins2, 1, replace=TRUE, prob=rep(1/n, n))

12. 


13. ## [1] 0

當(dāng)然最重要的一步是要能按照“時(shí)間”更新網(wǎng)絡(luò)節(jié)點(diǎn)被感染的信息。
?

1. keep =unlist(lapply(nearest_neighbors[,2], toss))

3. new_infected =as.numeric(as.character(nearest_neighbors[,1][keep >=1]))

5. diffusers =unique(c(as.numeric(diffusers), new_infected))

6. 


7. return(diffusers)}

8. 


9. set.seed(1);

?

開啟擴(kuò)散過程！

先看看S曲線吧：
?

1. # # "growth_curve"num_cum =unlist(lapply(1:i, function(x) length(infected[[x]]) ))

2. 


3. p_cum =num_cum time =1:i

4. 


5. ## Large initial population size (X=1000)

6. 


7. parms <-c(beta=0.01, gamma=0.1)

8. 


9. x0 <-c(S=49,I=1,R=0)a <-c("beta*S*I","gamma*I")

10. 


11. nu <-matrix(c(-1,0,+1,-1,0,+1),nrow=3,byrow=TRUE)

12. 


13. out <-ssa(x0,a,nu,parms,tf=4,simName="SIR model")

為了可視化這個(gè)擴(kuò)散的過程，我們用紅色來標(biāo)記被感染者。

1. # generate a palette#

2. 


3. plot(g, layout =layout.old)

4. 


5. set.seed(1)#

6. 


7. library(animation)# start the plot

8. 


9. m =1

1. same=numeric(0)

2. 


3. for(m in 2:length(health))

4. 


5. if(length(setdiff(health[[m ]],health[[m -1 ]]) )==0){same=c(same,m)

6. }

7. health=health[-same]

8. 


9. infected=infected[-same]#

如同在Netlogo里一樣，我們可以把網(wǎng)絡(luò)擴(kuò)散與增長(zhǎng)曲線同時(shí)展示出來：

1. set.seed(1)

2. 


3. # start the plot

4. 


5. m =1

6. 


7. p_cum=numeric(0)

8. 


9. h_cum=numeric(0)

10. 


11. i_cum=numeric(0)

12. 


13. while( m<50 ) {# start the plot

14. 


15. layout(matrix(c(1, 2, 1, 3), 2,2, byrow =TRUE), widths=c(3,1), heights=c(1, 1))

16. 


17. V(g)$color = "white"

18. 


19. V(g)$color[V(g)%in%infected[[m ]] ] = "red"

20. 


21. V(g)$color[V(g)%in%health[[m ]]] = "green"

22. 


23. if(m<=length(infected))

24. 


25. 


26. 


27. 


28. 


29. 


30. 


31. 


32. 


33. 


34. 


35. plot(pp~time, type ="h", ylab ="PDF", xlab ="Time",xlim =c(0,i), ylim =c(0,1), frame.plot =FALSE)

36. 


37. m =m +1

38. 


39. }

參考文獻(xiàn)

1.R語(yǔ)言泊松Poisson回歸模型分析案例

2.R語(yǔ)言進(jìn)行數(shù)值模擬：模擬泊松回歸模型

3.r語(yǔ)言泊松回歸分析

4.R語(yǔ)言對(duì)布豐投針（蒲豐投針）實(shí)驗(yàn)進(jìn)行模擬和動(dòng)態(tài)可視化

5.用R語(yǔ)言模擬混合制排隊(duì)隨機(jī)服務(wù)排隊(duì)系統(tǒng)

6.GARCH（1,1），MA以及歷史模擬法的VaR比較

7.R語(yǔ)言做復(fù)雜金融產(chǎn)品的幾何布朗運(yùn)動(dòng)的模擬

8.R語(yǔ)言進(jìn)行數(shù)值模擬：模擬泊松回歸模型

9.R語(yǔ)言對(duì)巨災(zāi)風(fēng)險(xiǎn)下的再保險(xiǎn)合同定價(jià)研究案例：廣義線性模型和帕累托分布Pareto distributions