# R/GWishart_NOij_Gibbs.R In NonDecompGraph: NonDecomposable Graph

#### Documented in GWishart_NOij_Gibbs

```#   Start from C\(i,j)
#
#     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
#     C: initial partial covariance matrix;
#     burnin, nmc : number of MCMC burnins and saved samples

GWishart_NOij_Gibbs <- function( bG, DG, adj, C, i, j, edgeij, burnin, nmc )
{
p <- dim( DG )[1]
## If no edge
if( edgeij == 0 ){
C[i,j]  <- 0
C[j,i]  <- 0
l       <- solve( DG[j,j,drop=FALSE] )
A       <- rwish( bG, l )
C_12    <- C[ j,-j,drop=FALSE]
C_22    <- C[-j,-j,drop=FALSE]
invC_22 <- solve( C_22 )
C[j,j]  <- A + C_12 %*% invC_22 %*% t( C_12 )
} else {
reorder   <- 1:p
reorder   <- c(reorder[-c( i, j )] ,i,j)
C_reorder <- C[reorder,reorder,drop=FALSE]
R         <- chol(C_reorder)
b         <- R[( p - 1 ),( p - 1 )]
m_post    <- -b * DG[i,j]/DG[j,j]
sig_post  <- 1/sqrt(DG[j,j])
R[p-1,p]  <- rnorm(1) * sig_post + m_post
R[p,p]    <- sqrt(rgamma(1, shape=bG/2, scale = 2/DG[j,j,drop=FALSE]) )  #gamrnd(bG/2,2/DG[j,j])
dR        <- dim(R)
C_updated <- t( R[,(dR[2]-1):dR[2],drop=FALSE]) %*% R[,dR[2],drop=FALSE]
C[i,j]    <- C_updated[1]
C[j,i]    <- C_updated[1]
C[j,j]    <- C_updated[2]
}
Sig = solve(C)

IsolatedNodeId <- which( colSums( adj ) == 1 ) # isolated node, i.e. nodes that are not connected to any other nodes
INsize         <- length( IsolatedNodeId )
if( ( burnin + nmc ) > 0 ){
for( iter in 1:( burnin + nmc ) ) {
###  Sample isolated nodes
if( INsize > 0 ){
for( i in 1:INsize ){
cliqueid <- IsolatedNodeId[i]
K_c      <- rwish( bG,solve( DG[cliqueid,cliqueid,drop=FALSE] ) )
C[cliqueid,cliqueid] <- K_c
Sig[cliqueid,cliqueid] <- 1/K_c
}
}
for( i in 1:( p - 1) ){
for( j in ( i + 1 ):p ){
if( adj[i,j] == 1 ){
cliqueid   <- c(i,j)
cliquesize <- 2
l          <- solve( DG[cliqueid,cliqueid,drop=FALSE] )
l          <- ( l + t( l ) )/2
A          <-  rwish(bG+cliquesize-1,l)
C_12       <-   C[ cliqueid,-cliqueid,drop=FALSE]
Sig_12     <- Sig[ cliqueid,-cliqueid,drop=FALSE]
Sig_22     <- Sig[-cliqueid,-cliqueid,drop=FALSE]
invSig_11  <- solve( Sig[cliqueid,cliqueid,drop=FALSE] )
invSig_11  <- (invSig_11 + t(invSig_11))/2
invC_22    <- Sig_22 - t( Sig_12 ) %*% invSig_11 %*% Sig_12
K_c        <- A + C_12 %*% invC_22 %*% t( C_12 )
K_c        <- ( K_c + t( K_c ) ) / 2 # Numerical stable
Delta      <- solve( C[cliqueid,cliqueid,drop=FALSE] - K_c )
C[cliqueid,cliqueid] <- K_c

## Rank 2 update Sig
Sig_bb <- Sig[cliqueid,cliqueid,drop=FALSE]
aa     <- solve( Delta - Sig_bb )
aa     <- ( aa + t( aa ) ) / 2
Sig    <- Sig + Sig[,cliqueid,drop=FALSE] %*% aa %*% t(Sig[,cliqueid,drop=FALSE])
}
}
}
}
}
return( list( C = C, Sig = Sig ) )
}
```

## Try the NonDecompGraph package in your browser

Any scripts or data that you put into this service are public.

NonDecompGraph documentation built on May 31, 2017, 4:56 a.m.