GWishart_BIPS_maximumClique: Sample from Gwishart distribution

In NonDecompGraph: NonDecomposable Graph

Description

Sample from Gwishart distribution

Usage

GWishart_BIPS_maximumClique(bG, DG, adj, C, burnin, nmc)

Arguments

bG	d.f.
DG	location
adj	adjacency matrix
C	initial partial covariance matrix
burnin	number of MCMC burnins
nmc	number of MCMC saved samples

Details

Sample C from Gwishart distribution with density: p(C) \propto |C|^{(bG-2)/2} exp(-1/2 tr(C DG)) where (1) bG : d.f. (2) DG: location (3) adj: adjacency matrix C: initial partial covariance matrix;

Value

C	Samples from G-Wishart distribution. The result is 3D array with c(dim(C)[1],dim(C)[2], nmc ).
Sig	Inverse of C The result is 3D array with c(dim(C)[1],dim(C)[2], nmc ).

Author(s)

Hao Wang ; Sophia Zhengzi Li

References

Wang and Li (2011) "Efficient Gaussian Graphical Model Determination without Approximating Normalizing Constant of the G-Wishart distribution " http://www.stat.sc.edu/~wang345/RESEARCH/GWishart/GWishart.html

Examples

  
##########################################################
##### e.g. (2) dense graph      ###########################
##########################################################
p     <- 10; 
alpha <- 0.3; 
J0    <- 0.5*diag(p); 
B0    <- 1*(matrix(runif(p*p),p,p)<alpha);
for( i in 1:p ){ 
    for( j in 1:i ){ 
        if ( B0[i,j] && i!=j ){
            J0[i,j] <- 0.5
        }
    }
}; 
J0    <- J0 + t(J0);
tmp   <- eigen( J0 ); w <- tmp$values;  
delta <- ( p * min( w ) - max( w ) ) / ( 1 - p ); 
J0    <- J0 + delta * diag( p );
D     <- diag(p) + 100*solve(J0);
DG    <- D;
adj   <- 1*(abs(J0)>1e-4); 
bG    <- 103;
    
burnin = 100; nmc = 100;

###  Maximum clique algorithm  ################
C = diag(p); # Initial value
resBIPSmax = GWishart_BIPS_maximumClique(bG,DG,adj,C,burnin,nmc);
C_save_maxC  <-resBIPSmax[[1]]
Sig_save_maxC<-resBIPSmax[[2]]

NonDecompGraph documentation built on May 2, 2019, 6:47 p.m.