原文連結：http://tecdat.cn/?p=14593
原文出處：拓端資料部落公衆号
與普通的擴散研究不同，網絡擴散開始考慮網絡結構對于擴散過程的影響。 這裡介紹一個使用R模拟網絡擴散的例子。
基本的算法非常簡單： 生成一個網絡:g(V, E)。 随機選擇一個或幾個節點作為種子（seeds）。 每個感染者以機率p（可視作該節點的傳染能力,通常表示為ββ）影響與其相連的節點。 其實這是一個最簡單的SI模型在網絡中的實作。S表示可感染（susceptible）, I表示被感染（infected）。易感态-感染态-恢複态(SIR)模型用以描述水痘和麻疹這類患者能完全康複并獲得終身免疫力的流行病。對于SIR流行病傳播模型，任意時刻節點隻能處于易感态(S)或感染态(I)或恢複态(R)。易感态節點表示未被流行病感染的個體，且可能被感染；感染态節點表示已經被流行病感染且具有傳播能力；恢複态節點則表示曾感染流行病且完全康複。與SIS模型類似，每一時間步内，每個感染态節點以機率λλ嘗試感染它的鄰居易感态節點，并以機率γγ變為恢複态。SIR模型可以表達為：
S = S（t）是易感個體的數量， I = I（t）是被感染的個體的數目， R = R（t）是恢複的個體的數目。
第二組因變量代表在三個類别的總人口的比例。是以，如果N是總人口（790萬在我們的例子），我們有
S（T）= S（T）/ N，人口的易感部分， Ⅰ（T）= I（t）的/ N的人口感染分數并 R（T）= R（t）的/ N，人口的康複部分。
解這個微分方程，我們可以得到累計增長曲線的表達式。有趣的是，這是一個logistic增長，具有明顯的S型曲線（S-shaped curve）特征。該模型在初期跨越臨界點之後增長較快，後期則變得緩慢。 因而可以用來描述和拟合創新擴散過程（diffusion of innovations）。 當然，對疾病傳播而言，SI模型是非常初級的（naive），主要因為受感染的個體以一定的機率恢複健康，或者繼續進入可以被感染狀态(S，據此擴充為SIS模型)或者轉為免疫狀态（R,據此擴充為SIR模型）。 免疫表示為R，用γγ代表免疫機率（removal or recovery rate)。對于資訊擴散而言，這種考慮暫時是不需要的。
第一步，生成網絡。
規則網
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)
g =connect.neighborhood(graph.ring(size), 2); plot(g) # 最近鄰耦合網絡
# 随機網絡g =erdos.renyi.game(size, 0.1)# 小世界網絡
g = rewire.edges(erdos.renyi.game(size, 0.1), prob = 0.8 )# 無标度網絡
g =barabasi.game(size) ; plot(g)
第二步，随機選取一個或n個随機種子。
# initiate the diffusers
seeds_num =1 diffusers =sample(V(g),seeds_num) ;
diffusers
## + 1/50 vertex:
## [1] 43
infected =list()
infected[[1]]=diffusers#
第三步，傳染能力
在這個簡單的例子中，每個節點的傳染能力是0.5，即與其相連的節點以0.5的機率被其感染,每個節點的回複能力是0.5，即其以0.5的機率被其回複。在R中的實作是通過抛硬币的方式來實作的。
## [1] 0
顯然，這很容易擴充到更一般的情況，比如節點的平均感染能力是0.128，那麼可以這麼寫： 節點的平均回複能力是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
當然最重要的一步是要能按照“時間”更新網絡節點被感染的資訊。
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);
開啟擴散過程！
先看看S曲線吧：
# # "growth_curve"num_cum =unlist(lapply(1:i, function(x) length(infected[[x]]) ))
p_cum =num_cum time =1:i
## Large initial population size (X=1000)
parms <-c(beta=0.01, gamma=0.1)
x0 <-c(S=49,I=1,R=0)a <-c("beta*S*I","gamma*I")
nu <-matrix(c(-1,0,+1,-1,0,+1),nrow=3,byrow=TRUE)
out <-ssa(x0,a,nu,parms,tf=4,simName="SIR model")
為了可視化這個擴散的過程，我們用紅色來标記被感染者。
# generate a palette#
plot(g, layout =layout.old)
set.seed(1)#
library(animation)# start the plot
m =1
same=numeric(0)
for(m in 2:length(health))
if(length(setdiff(health[[m ]],health[[m -1 ]]) )==0){same=c(same,m)
}
health=health[-same]
infected=infected[-same]#
如同在Netlogo裡一樣，我們可以把網絡擴散與增長曲線同時展示出來：
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
}
參考文獻
1.R語言泊松Poisson回歸模型分析案例
2.R語言進行數值模拟：模拟泊松回歸模型
3.r語言泊松回歸分析
4.R語言對布豐投針（蒲豐投針）實驗進行模拟和動态可視化
5.用R語言模拟混合制排隊随機服務排隊系統
6.GARCH（1,1），MA以及曆史模拟法的VaR比較
7.R語言做複雜金融産品的幾何布朗運動的模拟
8.R語言進行數值模拟：模拟泊松回歸模型
9.R語言對巨災風險下的再保險合同定價研究案例：廣義線性模型和帕累托分布Pareto distributions