library(showtext)
## Warning: package 'showtext' was built under R version 4.1.2
## Loading required package: sysfonts
## Loading required package: showtextdb
showtext_auto()
library(sand)
library(igraph)
set.seed(42)
g.er<-sample_gnp(100,0.02)
par(mfrow=c(1,1),mai=c(0.2,0.2,0.2,0.2))
plot(g.er,layout=layout_in_circle,vertex.size=3, vertex.label=NA,
main="E-R随机图:G(100,0.02)")
is_connected(g.er)
## [1] FALSE
#探讨G(100,0.02)的联通分支个数与各个联通分支节点规模
table(sapply(decompose(g.er), vcount))
##
## 1 2 3 4 71
## 15 2 2 1 1
#探索度分布、平均距离、直径以及聚集系数
mean(degree(g.er))
## [1] 1.9
hist(degree(g.er),col="lightblue",
xlab="Degree",ylab="Frequency",main="")
mean_distance(g.er)
## [1] 5.276511
diameter(g.er)
## [1] 14
transitivity(g.er)
## [1] 0.01639344
#探讨G(100,0.02)最大连通分支的平均距离、直径以及聚集系数
er.gc<-decompose(g.er)[[1]]
vcount(er.gc)
## [1] 71
mean_distance(er.gc)
## [1] 5.298592
diameter(er.gc)
## [1] 14
transitivity(er.gc) #计算聚集系数
## [1] 0.01685393
#计算G(100,0.02)最大连通分支边和节点联通度
edge_connectivity(er.gc)
## [1] 1
vertex_connectivity(er.gc)
## [1] 1
#探寻G(100,0.02)最大连通分支的点割集和边割集
er.gc.cut.vertex<-articulation_points(er.gc)
length(er.gc.cut.vertex)
## [1] 24
#计算G(100,0.02)最大连通分支laplacian矩阵特征值和特征向量
er.laplacian<-laplacian_matrix(g.er)
eig.anal <- eigen(er.laplacian)
eig.anal$values
## [1] 7.699012e+00 6.841590e+00 6.608220e+00 6.516514e+00 6.316720e+00
## [6] 6.165734e+00 5.484473e+00 5.264639e+00 5.009344e+00 4.947595e+00
## [11] 4.810788e+00 4.647818e+00 4.445484e+00 4.336075e+00 4.265182e+00
## [16] 4.127680e+00 4.000000e+00 3.913989e+00 3.862117e+00 3.643187e+00
## [21] 3.551151e+00 3.413502e+00 3.407506e+00 3.299268e+00 3.183980e+00
## [26] 3.024431e+00 3.000000e+00 3.000000e+00 2.862772e+00 2.810024e+00
## [31] 2.779206e+00 2.716778e+00 2.598796e+00 2.384501e+00 2.321351e+00
## [36] 2.298484e+00 2.130383e+00 2.000000e+00 2.000000e+00 1.866544e+00
## [41] 1.777191e+00 1.673121e+00 1.608747e+00 1.585163e+00 1.572716e+00
## [46] 1.444605e+00 1.373710e+00 1.328985e+00 1.319324e+00 1.104570e+00
## [51] 1.048841e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
## [56] 1.000000e+00 1.000000e+00 9.492613e-01 9.386918e-01 8.678446e-01
## [61] 7.521392e-01 7.312295e-01 6.961720e-01 6.701281e-01 6.255446e-01
## [66] 5.909537e-01 5.435254e-01 4.840486e-01 4.464249e-01 3.998715e-01
## [71] 3.817887e-01 3.647104e-01 3.216094e-01 2.822019e-01 2.200735e-01
## [76] 1.498285e-01 7.335784e-02 7.298415e-02 4.579743e-02 2.485383e-15
## [81] 1.842304e-15 1.229698e-15 5.772842e-16 4.189637e-16 3.397275e-16
## [86] 2.667933e-16 8.607625e-17 6.350102e-17 7.886900e-32 -5.421710e-32
## [91] -4.639810e-17 -7.188191e-17 -1.637232e-16 -1.832631e-16 -2.632191e-16
## [96] -3.812576e-16 -4.855212e-16 -8.202819e-16 -1.284250e-15 -2.833669e-15
par(mfrow=c(1,1))
plot(eig.anal$values, col="blue",
ylab="Eigenvalues of Graph Laplacian")
给定度序列,构建随机图。
degs<-c(2,2,2,2,3,3,3,3)
g1 <- sample_degseq(degs, method="vl")
g2 <- sample_degseq(degs, method="vl")
par(mfrow=c(1,2),mai=c(0.2,0.2,0.2,0.2))
plot(g1,vertex.label=NA)
plot(g2,vertex.label=NA)
#判断g1,g2是否同构
isomorphic(g1, g2)
## [1] FALSE
#计算图g1,g2的边数
c(ecount(g1),ecount(g2))
## [1] 10 10
生成一个和酵母菌蛋白质作用网络度序列相同的网络:
data(yeast)
degs<-degree(yeast)
fake.yeast<-sample_degseq(degs,method = "vl")
par(mfrow=c(1,2),mai=c(0.2,0.2,0.2,0.2))
plot(yeast,vertex.size=2,vertex.label=NA,main="yeast")
plot(fake.yeast,vertex.size=2,vertex.label=NA,main="fake.yeast")
#判断生成的网络和原来的酵母菌蛋白质相互作用网络度序列是否相同
all(degree(yeast)==degree(fake.yeast))
## [1] TRUE
#计算生成网络的聚集系数,并与原来的酵母菌蛋白质相互作用网络聚集系数比较
transitivity(yeast)
## [1] 0.4686178
transitivity(fake.yeast)
## [1] 0.04026804
is.connected(fake.yeast)
## [1] TRUE
diameter(fake.yeast)
## [1] 8
mean_distance(fake.yeast)
## [1] 3.494976
vertex_connectivity(fake.yeast)
## [1] 1
edge_connectivity(fake.yeast)
## [1] 1
上述结果表明新生成的网络大量聚类消失了
从排成一圈的\(N_{v}\)个节点集合开始,将每个节点与其中\(r\)个另邻居相连,之后,以 概率\(p\)移除每条边的一个端点,与另一个等概率选取的新节点相连,并避免生成自环和多重边。
#生成N=100,r=5,p=0.05的小世界网络
g.ws <- sample_smallworld(1, 100, 5, 0.05)
par(mfrow=c(1,2),mai=c(0.2,0.2,0.2,0.2))
plot(g.ws,vertex.size=2,layout=layout.circle,
vertex.label=NA,main="小世界网络")
#生成N=100,r=5,p=0网络
g.lat100 <- sample_smallworld(1, 100, 5, 0)
plot(g.lat100,vertex.size=2,layout=layout.circle,
vertex.label=NA,main="lattice 100")
#比较二者的聚集系数、直径、平均距离
c(transitivity(g.ws),transitivity(g.lat100))
## [1] 0.5051728 0.6666667
c(diameter(g.ws),diameter(g.lat100))
## [1] 5 10
c(mean_distance(g.ws),mean_distance(g.lat100))
## [1] 2.714545 5.454545
生成N=1000,r=10,p在很大范围内改变的小世界网络,说明即便p值很小,随机重连少量 边对于减小节点间距离的效果非常明显,同时还保留了相似的高聚集性:
steps<-seq(-4,-0.5,0.1)
len<-length(steps)
cl<-numeric(len)
md<-numeric(len)
ntrails<-100
for(i in (1:len)){
cl_trail<-numeric(ntrails)
md_trail<-numeric(ntrails)
for(j in (1:ntrails)){
g<-sample_smallworld(1,1000,10,10^steps[i])
cl_trail[j]<-transitivity(g)
md_trail[j]<-mean_distance(g)
}
cl[i]<-mean(cl_trail)
md[i]<-mean(md_trail)
}
plot(steps, cl/max(cl), ylim=c(0, 1), lwd=2, type="l",
col="blue", xlab=expression(log[10](p)),
ylab="Clustering and Average Path Length")
lines(steps, md/max(md),type="l", lwd=3, col="red")
模拟节点数N=100,每个新节点增加m=1条边的BA模型:
set.seed(42)
g.ba<-sample_pa(100,directed=FALSE)
par(mfrow=c(1,2))
plot(g.ba,vertex.size=2,layout=layout.circle,
vertex.label=NA,main="BA(100)")
hist(degree(g.ba),col="lightblue",xlab="Degree",
ylab="Frequency",main="Degree Distribution")
summary(degree(g.ba))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 1.00 1.00 1.98 2.00 9.00
par(mfrow=c(1,1),mai=c(0.2,0.2,0.2,0.2))
degree(g.ba)
## [1] 6 9 3 7 3 6 4 4 2 7 5 3 3 2 2 4 1 2 2 2 3 4 1 3 5 4 1 2 2 2 1 1 1 2 2 2 3
## [38] 2 2 2 1 3 1 4 2 1 2 2 2 1 1 3 1 1 1 1 1 2 1 1 2 2 1 1 1 1 1 1 1 1 1 1 1 1
## [75] 1 1 1 1 1 1 1 1 1 1 1 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1
dd.ba<-degree_distribution(g.ba)
dd.ba
## [1] 0.00 0.55 0.22 0.10 0.06 0.02 0.02 0.02 0.00 0.01
#计算平均距离、直径、聚集系数
mean_distance(g.ba)
## [1] 5.815556
diameter(g.ba)
## [1] 12
transitivity(g.ba)
## [1] 0
对空手道网络karate定义两个参考框架:
1、生成与空手道网络具有相同阶数(节点数)\(N_{v}=34\)和规模(边数)\(N_{e}=78\)的随机图;
2、生成与空手道网络具有相同度分布的随机图
重复实验1000次,对每个图使用相同的社团发现算法确定社团数量。
data(karate)
nv<-vcount(karate)
ne<-ecount(karate)
degs<-degree(karate)
ntrails<-1000
#生成具有和karate网络相同阶数和尺寸的经典随机图
num.comm.rg<-numeric(ntrails)
for(i in (1:ntrails)){
g.rg<-sample_gnm(nv,ne)
c.rg<-cluster_fast_greedy(g.rg)
num.comm.rg[i]<-length(c.rg)
}
#生成具有和karate网络相同度序列的广义随机图
num.comm.grg<-numeric(ntrails)
for(i in (1:ntrails)){
g.grg<-sample_degseq(degs,method="vl")
c.grg<-cluster_fast_greedy(g.grg)
num.comm.grg[i]<-length(c.grg)
}
#使用并列的条形图对结果进行对比
result<-c(num.comm.rg,num.comm.grg)
head(result)
## [1] 4 5 5 5 4 4
index<-c(rep(0,ntrails),rep(1,ntrails))
counts<-table(index,result)/ntrails
counts
## result
## index 3 4 5 6 7 8
## 0 0.006 0.256 0.517 0.198 0.022 0.001
## 1 0.011 0.227 0.564 0.189 0.009 0.000
barplot(counts, beside=TRUE, col=c("blue", "red"),
xlab="Number of Communities",
ylab="Relative Frequency",
legend=c("Fixed Size", "Fixed Degree Sequence"))
## 评估小世界性 igraphdata中的macaque数据集包含了一个猕猴大脑区域视觉皮层触觉功能网络。
定义有向图聚集系数: \[cl(v)=\frac{(A+A^{T})_{vv}^{3}}{2[d_{v}^{tot}(d_{v}^{tot}-1)-2(A^2)_{vv}]}\] 其中A为邻接矩阵,\(d_{v}^{tot}\)是节点\(v\)的总度值(入度加出度)。 可以自编有向图聚集系数函数如下:
clust.coef.dir<-function(graph){
A<-as.matrix(get.adjacency(graph))
S<-A+t(A)
deg<-degree(graph,mode=c("total"))
num<-diag(S %*% S %*% S) #分子
denom<-diag(A %*% A)
denom<-2*(deg*(deg-1)-2*denom)
cl<-mean(num/denom)
return(cl)
}
模拟抽取与macaque网络同阶和同规模的有向经典图并评估每个图的聚集系数和平均有向路径长度,并与 原来的macaque网络进行对比分析,实验重复1000次:
library(igraphdata)
data(macaque)
summary(macaque)
## IGRAPH f7130f3 DN-- 45 463 --
## + attr: Citation (g/c), Author (g/c), shape (v/c), name (v/c)
nv<-vcount(macaque)
ne<-ecount(macaque)
plot(macaque,vertex.size=3,vertex.label=NA)
ntrails<1000
## [1] FALSE
cl.rg<-numeric(ntrails)
md.rg<-numeric(ntrails)
for(i in (1:ntrails)){
g.rg<-sample_gnm(nv,ne,directed = TRUE)
cl.rg[i]<-clust.coef.dir(g.rg)
md.rg<-mean_distance(g.rg)
}
summary(cl.rg)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.2145 0.2301 0.2339 0.2339 0.2377 0.2517
summary(md.rg)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.836 1.836 1.836 1.836 1.836 1.836
c(clust.coef.dir(macaque),mean(cl.rg))
## [1] 0.5501073 0.2339433
c(mean_distance(macaque),mean(md.rg))
## [1] 2.148485 1.835859