This function estimates the diagonal of the precision matrix by symmetry-enforced likelihood minimization using shortest path trees, when the true value of the coefficient matrix
B is known or has already been estimated. The observations of the data matrix
X are assumed to have zero mean.
The data matrix.
The coefficient matrix.
The method to select the shortest path trees to be considered.
The matrix orresponding to outliers.
The shortest path trees for each connected component could be build under one of the following specifications :
only the presence or absence of edges is considered to build the shortest path trees,
the root of each shortest path tree is chosen as the node of maximal degree among those in the connected component,
the root of each shortest path tree - of maximum height equal to 1 - is chosen as the node of maximal degree,
for each connected component, the selected shortest path tree is the maximum weighted tree among all shortest path trees.
When Theta is not NULL, we consider an additive contamination model. We assume that X = Y + E is observed, denoting the outlier-free data by Y and the matrix of errors by E. In this case, the matrix Theta should be equal to E * B.
This function returns the diagonal of the precision matrix associated with X as a vector.
Arnak Dalalyan and Samuel Balmand.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
## build the true precision matrix p <- 12 # number of variables Omega <- 2*diag(p) Omega[1,1] <- p Omega[1,2:p] <- 2/sqrt(2) Omega[2:p,1] <- 2/sqrt(2) ## compute the true diagonal of the precision matrix Phi <- 1/diag(Omega) ## generate the design matrix from a zero-mean Gaussian distribution require(MASS) n <- 100 # sample size X <- mvrnorm(n,rep.int(0,p),ginv(Omega)) ## compute the sample mean barX <- apply(X,2,mean) ## subtract the mean from all the rows X <- t(t(X)-barX) ## estimate the coefficient matrix B <- DESP_SRL_B(X,lambda=sqrt(2*log(p))) ## compute the squared partial correlations SPC <- DESP_SqPartCorr(B,n) ## re-estimate the coefficient matrix by ordinary least squares B_OLS <- DESP_OLS_B(X,SPC) ## estimate the diagonal of the precision matrix and get its inverse hatPhiSPT <- 1/DESP_SPT(X,B_OLS,method=2.1) ## measure the performance of the estimation using l2 vector norm sqrt(sum((Phi-hatPhiSPT)^2))