Reputation: 11
I am working on R for generating graph energy of rank-1 heterogeneous random graph in which edges sharing same nodes have 0.5 correlation and edges with different nodes are mutually exclusive. EN1 is the energy of graph but gives all 0 values. Where is the mistake?
library(MASS)
#library(EnvStats)
#library(energy)
library(mvtnorm)
n <- 3
Pn <- 0.1
kappa <- 1
r <- 0.5
M <- 20
EN1 <- numeric(M)
for (rep in 1:M) {
A <- matrix(0, n, n)
for (i in 1:(n - 1)) {
for (j in (i + 1):n)
for (k in (j + 1):n)
for (l in (k + 1):n)
A[i, j] <- rbinom(1, 1, Pn*exp(-kappa*i/n)*exp(-kappa*j/n))
A[k, l] <- rbinom(1, 1, Pn*exp(-kappa*k/n)*exp(-kappa*l/n))
}
A <- A + t(A)
A
}
for (rep in 1:M) {
cov_matrix <- matrix(r, nrow=n, ncol=n)
diag(cov_matrix) <- 0
for (i in 1:(n - 1)) {
for (j in (i + 1):n)
for (k in (j + 1):n)
for (l in (k + 1):n)
B <- mvtnorm::rmvnorm(n, sigma=cov_matrix)
B[i, j] <- rbinom(1, 1, Pn*exp(-kappa*i/n)*exp(-kappa*j/n))
B[k, l] <- rbinom(1, 1, Pn*exp(-kappa*k/n)*exp(-kappa*l/n))
}
B <- B + t(B)
B
}
C <- A + B
C <- C + t(C)
C
ei1 <- eigen(C, symmetric=TRUE, only.values=TRUE)
ei1$values
EN1[rep] <- sum(abs(ei1$values))
EN1
mean(EN1) # report mean
sd(EN1) # report standard deviation
Upvotes: 0
Views: 42
Reputation: 19339
It looks like you want to put this code block:
C <- A + B
C <- C + t(C)
ei1 <- eigen(C, symmetric=TRUE, only.values=TRUE)
EN1[rep] <- sum(abs(ei1$values))
within the rep
for loop, like so:
for (rep in 1:M) {
cov_matrix <- matrix(r, nrow=n, ncol=n)
diag(cov_matrix) <- 0
for (i in 1:(n - 1)) {
for (j in (i + 1):n)
for (k in (j + 1):n)
for (l in (k + 1):n)
B <- mvtnorm::rmvnorm(n, sigma=cov_matrix)
B[i, j] <- rbinom(1, 1, Pn*exp(-kappa*i/n)*exp(-kappa*j/n))
B[k, l] <- rbinom(1, 1, Pn*exp(-kappa*k/n)*exp(-kappa*l/n))
}
B <- B + t(B)
C <- A + B
C <- C + t(C)
ei1 <- eigen(C, symmetric=TRUE, only.values=TRUE)
EN1[rep] <- sum(abs(ei1$values))
}
This gives:
> EN1
[1] 7.076329 1.880778 7.902730 2.055223 4.731054 6.326656 10.886720 7.226352
[9] 7.703960 8.511892 6.168880 4.670598 3.709914 5.671895 1.671177 15.059931
[17] 11.200192 4.440183 3.822154 3.661713
> mean(EN1) # report mean
[1] 6.218917
> sd(EN1) # report standard deviation
[1] 3.405246
Upvotes: 0