library(igraph)
library(Biodem)
library(ggplot2)
library(gridExtra)
b=c(0,0.9,0.1,0.1,0.9,0.9,0,0.9,0.1,0.9,0.1,0.9,0,0.1,0.9,0.1,0.1,0.1,0,0.9,0.9,0.9,0.9,0.9,0)
a=matrix(b,ncol = 5,byrow = T)
a[lower.tri(a)]=0#将下三角置为0
d=data.frame(stringsAsFactors = F)
for(i in 1:5){
p1=rep(i,5)
p2=c(1:5)
weight=a[i,]
temp=data.frame(p1,p2,weight)
d=rbind(d,temp)
}
d1=d[d[,3]!=0,]
g = graph.data.frame(d1, directed = F)
plot(g,edge.width=d1$weight*3,edge.color='black',layout=layout.circle(g))
fun_WUN=function(a){
n=ncol(a)
one=matrix(rep(1,n),ncol = 1)
zero=one%*%t(one)-diag(n)
x=seq(0,1,0.001)
C_i_t=matrix(nrow = length(x),ncol = n)
for(j in 1:length(x)){
at=matrix(nrow=n,ncol=n)
at[a>x[j]]=1
at[is.na(at)]=0
gamma_t=mtx.exp(at,3)
Gamma_t=at%*%zero%*%at
for(i in 1:n){
C_i_t[j,i]=ifelse(Gamma_t[i,i]==0,0,gamma_t[i,i]/Gamma_t[i,i])
}
}
xx=data.frame(t=x,y = C_i_t,stringsAsFactors = F)
return(xx)
}
b=c(0,0.9,0.1,0.1,0.9,0.9,0,0.9,0.1,0.9,0.1,0.9,0,0.1,0.9,0.1,0.1,0.1,0,0.9,0.9,0.9,0.9,0.9,0)
a=matrix(b,ncol = 5,byrow = T)
xx=fun_WUN(a)
g=ggplot(xx)
g1=g+geom_line(mapping = aes(xx$t,xx$y.1))+xlab("t")+ylab("C1(t)")+ggtitle("Nodes 1 and 3")
g2=g+geom_line(mapping = aes(xx$t,xx$y.2))+xlab("t")+ylab("C2(t)")+ggtitle("Nodes 2")
g4=g+geom_line(mapping = aes(xx$t,xx$y.4))+xlab("t")+ylab("C4(t)")+ggtitle("Nodes 4")
g5=g+geom_line(mapping = aes(xx$t,xx$y.5))+xlab("t")+ylab("C5(t)")+ggtitle("Nodes 5")
grid.arrange(g1, g2, g4, g5, nrow=2, ncol=2)
C_i=(colSums(xx)/nrow(xx))[-1]
C=mean(C_i)
print(C_i)#clustering coefficient for node i
## y.1 y.2 y.3 y.4 y.5
## 0.8991009 0.6327006 0.8991009 0.0999001 0.3663004
print(C)#clustering coefficient
## [1] 0.5794206
下图是论文中的结果
结果跟论文不一样啊。。可是论文的结果明显不符合它的数据
fun_WDN8=function(a){
n=ncol(a)
one=matrix(rep(1,n),ncol = 1)
zero=one%*%t(one)-diag(n)
x=seq(0,1,0.001)
C_i_t=matrix(nrow = length(x),ncol = n)
for(j in 1:length(x)){
at=matrix(nrow=n,ncol=n)
at[a>x[j]]=1
at[is.na(at)]=0
gamma_t=mtx.exp(at+t(at),3)/8
Gamma_t=(at+t(at))%*%zero%*%(at+t(at))/4
for(i in 1:n){
C_i_t[j,i]=ifelse(Gamma_t[i,i]==0,0,gamma_t[i,i]/Gamma_t[i,i])
}
}
xx=data.frame(t=x,y = C_i_t,stringsAsFactors = F)
return(xx)
}
same as WUN
Evolution of Ci(t) as t ranges from 0 to 1 for one node in a WUN consisting of 100 nodes with edge weights drawn from a uniform distribution on (0, 1).
a=matrix(runif(10000,min=0,max=1), ncol = 100, nrow=100, byrow = T)
a[lower.tri(a)]=0
diag(a)=0
a=t(a)+a
xx=fun_WUN(a)
g=ggplot(xx)
g1=g+geom_line(mapping = aes(xx$t,xx$y.1))+xlab("t")+ylab("C1(t)")+ggtitle("Nodes 1")
print(g1)
跟论文还是不太一样- -
C_i=(colSums(xx)/nrow(xx))[-1]
C=mean(C_i)
print(C_i)#clustering coefficient for node i
## y.1 y.2 y.3 y.4 y.5 y.6 y.7
## 0.4959219 0.5010233 0.4874833 0.5090941 0.5052689 0.5037506 0.5119989
## y.8 y.9 y.10 y.11 y.12 y.13 y.14
## 0.5017656 0.5051738 0.5017891 0.5090628 0.5049541 0.4875044 0.5090246
## y.15 y.16 y.17 y.18 y.19 y.20 y.21
## 0.5259459 0.5188837 0.5141408 0.5100411 0.5053733 0.5062289 0.4970485
## y.22 y.23 y.24 y.25 y.26 y.27 y.28
## 0.4916233 0.5050157 0.5160245 0.4907722 0.4962520 0.5053985 0.5040177
## y.29 y.30 y.31 y.32 y.33 y.34 y.35
## 0.5034800 0.5205683 0.5153651 0.4951459 0.4956400 0.4863941 0.4935662
## y.36 y.37 y.38 y.39 y.40 y.41 y.42
## 0.5145860 0.4948101 0.4862238 0.5060069 0.4984667 0.4885635 0.5007670
## y.43 y.44 y.45 y.46 y.47 y.48 y.49
## 0.5078548 0.5097353 0.5178197 0.5050618 0.5007630 0.5146748 0.4976050
## y.50 y.51 y.52 y.53 y.54 y.55 y.56
## 0.4960986 0.4989776 0.5071091 0.4992430 0.4990020 0.4925128 0.5003875
## y.57 y.58 y.59 y.60 y.61 y.62 y.63
## 0.5105412 0.5010946 0.4992111 0.5086728 0.4981406 0.4999308 0.5150641
## y.64 y.65 y.66 y.67 y.68 y.69 y.70
## 0.4872839 0.5198870 0.5036227 0.5073775 0.5048351 0.4870332 0.5070581
## y.71 y.72 y.73 y.74 y.75 y.76 y.77
## 0.5062427 0.5187670 0.5163321 0.4782912 0.4992425 0.5086012 0.5008562
## y.78 y.79 y.80 y.81 y.82 y.83 y.84
## 0.5122595 0.4935235 0.5025948 0.5084875 0.5000987 0.4957769 0.4965781
## y.85 y.86 y.87 y.88 y.89 y.90 y.91
## 0.4895896 0.4893925 0.4930107 0.5203137 0.4995830 0.4980117 0.4914758
## y.92 y.93 y.94 y.95 y.96 y.97 y.98
## 0.5108642 0.5237322 0.5075171 0.5012531 0.5024007 0.5050857 0.5013012
## y.99 y.100
## 0.5031993 0.5003922
print(C)#clustering coefficient
## [1] 0.5029154
t=read.csv("42部门投入产出表.csv",header=T,stringsAsFactors = F)[,-1]
X=sapply(t,sum)
#求直接消耗系数a
a=matrix(nrow = 42,ncol = 42)
for(i in 1:42){
for(j in 1:42){
a[i,j]=t[i,j]/X[j]
}
}
#标准化a
w=a/max(a)
xx=fun_WDN8(w)
g=ggplot(xx)
g1=g+geom_line(mapping = aes(xx$t,xx$y.1))+xlab("t")+ylab("C1(t)")+ggtitle("Nodes 1")
g2=g+geom_line(mapping = aes(xx$t,xx$y.2))+xlab("t")+ylab("C2(t)")+ggtitle("Nodes 2")
g4=g+geom_line(mapping = aes(xx$t,xx$y.3))+xlab("t")+ylab("C3(t)")+ggtitle("Nodes 3")
g5=g+geom_line(mapping = aes(xx$t,xx$y.4))+xlab("t")+ylab("C4(t)")+ggtitle("Nodes 4")
grid.arrange(g1, g2, g4, g5, nrow=2, ncol=2)
C_i=(colSums(xx)/nrow(xx))[-1]
C=mean(C_i)
yy=data.frame(C_i,1:42,text4)
names(yy)=c("C_i","No.","Dep")
print(yy[sort.list(yy[,1]),])
## C_i No. Dep
## y.38 0.06660388 38 居民服务、修理和其他服务
## y.42 0.06695306 42 公共管理、社会保障和社会组织
## y.40 0.07052497 40 卫生和社会工作
## y.3 0.07116783 3 石油和天然气开采产品
## y.37 0.08477626 37 水利、环境和公共设施管理
## y.5 0.08663060 5 非金属矿和其他矿采选产品
## y.39 0.08732557 39 教育
## y.29 0.10344554 29 批发和零售
## y.33 0.10471835 33 金融
## y.27 0.10709412 27 水的生产和供应
## y.36 0.12027013 36 科学研究和技术服务
## y.22 0.12203521 22 其他制造产品
## y.35 0.12649177 35 租赁和商务服务
## y.24 0.13477441 24 金属制品、机械和设备修理服务
## y.41 0.13798905 41 文化、体育和娱乐
## y.10 0.14800081 10 造纸印刷和文教体育用品
## y.11 0.15309952 11 石油、炼焦产品和核燃料加工品
## y.18 0.16875060 18 交通运输设备
## y.34 0.17116720 34 房地产
## y.26 0.20420247 26 燃气生产和供应
## y.17 0.20770153 17 专用设备
## y.23 0.22714860 23 废品废料
## y.9 0.24023710 9 木材加工品和家具
## y.28 0.24080627 28 建筑
## y.15 0.26048516 15 金属制品
## y.30 0.26891728 30 交通运输、仓储和邮政
## y.19 0.27449010 19 电气机械和器材
## y.4 0.27666028 4 金属矿采选产品
## y.8 0.28322996 8 纺织服装鞋帽皮革羽绒及其制品
## y.32 0.28379768 32 信息传输、软件和信息技术服务
## y.31 0.31463373 31 住宿和餐饮
## y.25 0.32337877 25 电力、热力的生产和供应
## y.21 0.32745802 21 仪器仪表
## y.20 0.33761508 20 通信设备、计算机和其他电子设备
## y.16 0.41483758 16 通用设备
## y.6 0.41570818 6 食品和烟草
## y.13 0.43885788 13 非金属矿物制品
## y.2 0.46474023 2 煤炭采选产品
## y.14 0.46557297 14 金属冶炼和压延加工品
## y.1 0.48417152 1 农林牧渔产品和服务
## y.7 0.74261815 7 纺织品
## y.12 0.79784792 12 化学产品
print(C)#clustering coefficient
## [1] 0.2482604