require(graph)
require(mvtnorm)
nVar <- 50 ## number of variables
nObs <- 100 ## number of observations to simulate
set.seed(123)
g <- randomEGraph(as.character(1:nVar), p=0.15)
Sigma <- qpG2Sigma(g, rho=0.5)
X <- rmvnorm(nObs, sigma=as.matrix(Sigma))
## MLE of the sample covariance matrix
S <- cov(X)
## more efficient MLE of the sample covariance matrix using IPF
clqs <- qpGetCliques(g, verbose=FALSE)
S_ipf <- qpIPF(S, clqs)
## get the adjacency matrix and put the diagonal to one
A <- as(g, "matrix")
diag(A) <- 1
## entries in S and S_ipf for present edges in g should coincide
max(abs(S_ipf[A==1] - S[A==1]))
## entries in the inverse of S_ipf for missing edges in g should be zero
max(solve(S_ipf)[A==0])
Run the code above in your browser using DataLab