# R/log_GWishart_NOij_pdf.R In NonDecompGraph: NonDecomposable Graph

#### Documented in log_GWishart_NOij_pdf

```# Compute log p(C\c(i,j) ) upto the normalizing constant of G-Wishart
log_GWishart_NOij_pdf<-function(b_prior,D_prior,C,i,j,edgeij)
{
if( edgeij ==0 ){
C[i,j]    <- 0;
C[j,i]    <- 0;
C_12      <- C[j,-j,drop=FALSE];  #C_12[,j] = [];
C_22      <- C[-j,-j,drop=FALSE]; #C_22[j,] = []; C_22[,j] = [];
invC_22   <- solve(C_22);
C_new     <- C;
C_new[j,j]<- C_12%*%invC_22%*%t(C_12);
l         <- -log_iwishart_InvA_const( b_prior,D_prior[j,j]) + ( b_prior - 2 )/2 * log( det( C_22 ) ) -sum(diag( D_prior %*% C_new ))/2;
} else {

cliqueid  <- c(i,j);
C_12      <- C[cliqueid,-cliqueid,drop=FALSE];  #C_12[,cliqueid]=[];
C_22      <- C[-cliqueid,-cliqueid,drop=FALSE]; #C_22[cliqueid,]=[]; C_22[,cliqueid]=[];
invC_22   <- solve( C_22 );
A         <- C[cliqueid,cliqueid,drop=FALSE] - C_12 %*% invC_22 %*% t( C_12 );
A         <- ( t( A ) + A ) / 2;
log_Joint <- (b_prior-2)/2*log(det(C))-sum(diag(D_prior%*%C))/2;
logK2by2  <- log_GWishart_pdf(A,b_prior,D_prior[cliqueid,cliqueid,drop=FALSE],matrix(1,2,2),10,1);
V         <- solve( D_prior[cliqueid,cliqueid,drop=FALSE] );
D_priorii <- 1/V[1,1,drop=FALSE];
logKii    <- log_GWishart_pdf( A[1,1,drop=FALSE], b_prior + 1, D_priorii, matrix(1,1,1), 10, 1 );
l         <- log_Joint+logKii-logK2by2;

} #end
return(l)
}

# function l = log_GWishart_NOij_pdf(b_prior,D_prior,C,i,j,edgeij)

# % Compute log p(C\c(i,j) ) upto the normalizing constant of G-Wishart

# if edgeij ==0

#         C_12 = C(j,:);  C_12(:,j)=[];
#         C_22 = C; C_22(j,:)=[]; C_22(:,j)=[];
#         invC_22 = inv(C_22);

#         c = C_12*invC_22*C_12';

# %        A = 1;
# %        C(j,j) = A + C_12*invC_22*C_12';
# %
# %        log_Joint = (b_prior-2)/2*log(det(C))-trace(D_prior*C)/2;
# %        logKii = log_GWishart_pdf(A,b_prior,D_prior(j,j),1,1,1);
# %        l = log_Joint-logKii;

#        C_new = C; C_new(j,j)=c;
#        l = -log_iwishart_InvA_const(b_prior,D_prior(j,j))+(b_prior-2)/2*log(det(C_22))-trace(D_prior*C_new)/2;

# else

#         C_12 = C(cliqueid,:);  C_12(:,cliqueid)=[];
#         C_22 = C; C_22(cliqueid,:)=[]; C_22(:,cliqueid)=[];
#         invC_22 = inv(C_22);
#         A = C(cliqueid,cliqueid)-C_12*invC_22*C_12';        A = (A'+A)/2;

#         log_Joint = (b_prior-2)/2*log(det(C))-trace(D_prior*C)/2;

#         logK2by2 = log_GWishart_pdf(A,b_prior,D_prior(cliqueid,cliqueid),ones(2),10,1);

#         V = inv(D_prior(cliqueid,cliqueid)); D_priorii = 1/V(1,1);
#         logKii = log_GWishart_pdf(A(1,1),b_prior+1,D_priorii,1,10,1);

#         l = log_Joint+logKii-logK2by2;

# end
```

## 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.