ランダムグラフの構成単体数分布 - ryamadaの遺伝学・遺伝統計学メモ

こちらでDEC(Discrete Exterior Calculus)をいじっている
その冒頭で、グラフの隣接行列から、行列演算を繰り返すことで、グラフにあるp=1,2,3... 単体を列挙する話がある
Rで書いてみる
エッジ存在率を定めて発生させるランダムグラフがあったときに、単体の数がどのように分布するかをシミュレーションしてみる
今のランダムグラフ発生方法だと、1-単体(辺)は二項分布になる。p:(1-p)
2^単体(三角形)の場合は、p^3:(1-p^3)の二項分布…というように単体の構成辺数乗に関する二項分布になるはず(だけれども、比較的ノード数が少ない場合には、単体数の平均は二項分布予想とよく合うが、分散は、合わない。それは、「取りやすい値」に実単体数が集中するためのようだ)
やってみる

library(igraph)
n <- 50
p <- 0.5
n.iter <- 1000
res <- matrix(0,n.iter,n)
for(ii in 1:n.iter){
E <- matrix(0,n,n)
E[which(upper.tri(E))] <- sample(0:1,n*(n-1)/2,replace=TRUE,prob=c(1-p,p))
g <- graph.adjacency(E)
#plot(g)
#el <- get.edgelist(g)
S.list <- list()
S.list[[1]] <- diag(rep(1,n))
cnt <- 1
loop <- TRUE
while(loop){
	if(length(S.list[[cnt]][,1])==0){
		break
	}
	tmp <- S.list[[cnt]] %*% E
	tmp2 <- which(tmp==cnt,arr.ind=TRUE)
	if(length(tmp2[,1])>0){
		cnt <- cnt
		S.list[[cnt+1]] <- matrix(0,length(tmp2[,1]),n)
		for(i in 1:length(tmp2[,1])){
			u <- tmp2[i,1]
			v <- tmp2[i,2]
			pre <- which(S.list[[cnt]][u,]==1)
			S.list[[cnt+1]][i,c(pre,v)] <- 1
		}
	}else{
		loop <- FALSE
	}
	cnt <- cnt+1
}

tmp <- sapply(S.list,nrow)
res[ii,1:length(tmp)] <- tmp
}

boxplot(res)
# 生成されうる最大の単体数(対数)
s.full <- lchoose(n,1:n)
# 生成された単体数の平均と分散
s.mean <- apply(res,2,mean)
s.var <- apply(res,2,var)

plot(s.mean)
plot(s.var)
# 単体構成エッジ数
n.e <- (1:n)*((1:n)-1)/2
# 単体構成エッジ数を考慮した、単体の生起確率
p. <- p^n.e
# 二項分布からの理論的平均数と分散
t.mean <- exp(s.full)*p.
t.var <- exp(s.full)*p.*(1-p.)

plot(s.mean,t.mean)
abline(0,1,col=2)

plot(s.var,t.var)
abline(0,1,col=2)