bdmcmc: BDMCMC algorithm for undirected graph estimation

Description

This function is the BDMCMC algorithm for model selection in undirected Gaussian graphical models based on birth-death MCMC methodology. It works for both low and high dimensional graphs.

Usage

bdmcmc(data, n = NULL, meanzero = FALSE, iter = 5000, burn = floor(iter/2), 
       skip = 1, gamma.b = 1, prior.g = "Uniform", b = 3, D = NULL, 
	       A = "full", MCiter = 10, summary = FALSE, verbose = TRUE, all.A = FALSE)

Arguments

data

it could be a matrix or a data.frame of the data ( $n \times p$ matrix) or it could be a covariance matrix as $S=x'x$ which $x$ is a matrix of the data.

the number of observations.

meanzero

logical: if TRUE that means the data have zero mean; if FALSE (default) the mean of the data should not be zero.

iter

the number of iterations for the BDMCMC algorithm.

burn

the number of burn-in iterations for the BDMCMC algorithm.

skip

it is the option for regularly saving part of the iterations. Default is 1 (for saving all iterations).

gamma.b

the birth rates for the birth-death process.

prior.g

a character for selecting a prior distribution for the graph. It can be either "Uniform" (default) or "Poisson". "Uniform" means discrete uniform distribution for prior distribution of the graph ($G \sim DU(\mathcal{G})$ in which $\mathcal{G}$ is all p

the degree of freedom for G-Wishart distribution, $W_G(b,D)$. G-Wishart distribution is the prior distribution of precision matrix.

the positive definite matrix for G-Wishart distribution, $W_G(b,D)$. G-Wishart distribution is the prior distribution of precision matrix.

matrix A shows the starting graph for BDMCMC algorithm. It is an upper triangular matrix in which $a_{ij}=1$ if there is a link between nodes $i$ and $j$, otherwise $a_{ij}=0$. It can be either "full" (default) or "empty" or "glasso". "full" means we s

MCiter

the number of iterations for Monte Carlo approximation of normalizing constant in G-Wishart distribution. We need it only when number of nodes in the graph is less than 8.

summary

logical: if TRUE you will see the summary result of BDMCMC algorithm; if FALSE (default), the output will be a list as below.

verbose

logical: if verbose = TRUE, tracing number of iterations in BDMCMC algorithm printing is abled. The default value is TRUE.

all.A

logical: if all.A = FALSE, we save the matrix A after burn-in iteration. The default value is FALSE. If all.A = FALSE, we save all matrix A from starting point of iteration.

Value

Returns an output of the BDMCMC algorithm which is list-like and contains the following:
Asa vector which includes the adjacency matrix for all iterations after burn-in iterations.
lambdaa vector which includes the waiting times for all iterations after burn-in iterations.
allAa vector which includes all adjacency matrix for all iterations (includes burn-in iterations). We need the 'allA' for checking the convergency of the BDMCMC algorithm.

References

Mohammadi, A. and E. C. Wit (2012). Gaussian graphical model determination based on birth-death MCMC inference, arXiv:1210.5371v4. http://arxiv.org/abs/1210.5371v4

Examples

Run this code

p <- 8 # number of nodes 
  # "truK" is the precision matrix of true graph
  truK <- diag(p)
  for (i in 1:(p-1)) truK[i,i+1] <- truK[i+1,i] <- 0.5
  truK[1,p] <- truK[p,1] <- 0.4
  truK # precision matrix of the true graph
  
  # generate the data (200 observations) from multivariate normal 
  # distribution with mean zero and percision matrix "truK"
  data <- mvrnorm(200, c(rep(0,p)), solve(truK))
  
  bdmcmc(data, meanzero = T, iter = 2000, summary = T)