
\documentclass[12pt]{article}
\title{Notes for Different Graphical Gaussian Models}
\author{Guiyuan Lei\footnote{Guiyuan.Lei@ncl.ac.uk}}
\date{September, 2006}
\begin{document}
\maketitle
\section{HdBCS: Bayesian Covariance Selection in High-Dimensions}
The model of Bayesian Covariance Selection in High-Dimensions 
\cite{Dobra_2004} is designed to perform covariance selection for datasets with tens or possibly hundreds of thousands of variables. 
\subsection{The algorithm}
The main idea of HdBCS is to identify good regression models for each variable and to combine these models in a joint multivariate normal distribution using the chain rule (which is essentially equivalent to constructiong a DAG or Bayesian network). The algorithm starts with dependency network which is built from sets of regressions seperately, then generates an appropriate ordering of the variables to underlie a compositional network. With this ordering, variable selection for each linear model for the ordered variables produces a DAG. Finally, by removing arrows and marrying all parents of any node in this DAG, an undirected graph is delivered.\\

For analysis and selection of regressions, the model utilizes the inverse Wishart prior, applies for forward/backward varialbe selection with respect to the defined set of candidate predictors, and score each regression by the posterior model probabilities. See section 3 in the paper \cite{Dobra_2004} for detail. \\
%equation (4) is marginal model likelihood

Once the predictors of each variable are obtained, the order of all variables should be generated to construct compositional network. The strategy is to: the ordering aims to maximize the resulting posterior probability on the resulting model in an iterative process that ,at each step, decreases the overall score (equation 5 in the paper) the least. See section 4.1 for fully described construction.

\subsection{The software}
There is one software implemented for this model. Codes can be downloaded from 
http://www.stat.duke.edu/$\sim$adobra/hdbcs.html\\
%http://www.stat.duke.edu/~adobra/hdbcs.html
%HdBCS: Bayesian Covariance Selection in High-Dimensions
%Adrian Dobra (adobra@stat.duke.edu)

The search procedure has three separate steps: filtering, generating good starting models and improving the starting models until convergence. \\ 

The fisrt step is to select a relatively rich set of possible predictors for each variable. The search strategy is forward/backward procedure. This procedure is for each variable seperately. Starting with random selected predictors for each variable $x_i$, consider the edges between all other variables $x_j, j=[-i]$ with this variable, if $x_j$ is already predictor of variable $x_i$, propose to delete $x_j$ as predictor of $x_i$, otherwise, propose to add $x_j$ as its predictor. Calculate the ratio of posterior probability between proposed model and current model (one predictor difference between these two models), and accept/reject the proposed model according to the accept propobility. Repeat the above process for many iterations to get set of possible predictors for each variable. The pseudo code for the search procedure can be describled as following:\\

%see gibbsreg.cpp
For each variable $x_i$, 
\begin{enumerate}
\item Initial model, randomly select predictors of $x_i$ from all other variables $x_j, j=$\{$1,\cdots,p$\}$\setminus$\{$i$\} 
\item Iteration step: repeat the following process for a fixed times:
 \begin{enumerate}
 \item[] For all other variables $x_j, j=$\{$1,\cdots,p$\}$\setminus$\{$i$\}, randomly select $x_j$ (each variable can be selected once)
       \begin{enumerate}
        \item Propose new model from current model    
          \begin{enumerate}
             \item[] If ($x_j$ is not currently predictor of $x_i$), calculate the log posterior probability with $x_j$ is   predictor of $x_i$
             \item[] else, calculate the log posterior probability calculate the log posterior probability with $x_j$ is not predictor of $x_i$
           \end{enumerate}
         \item Calculate accept propobility and accept/reject proposed model, if $x_j$ is a predictor in new model, update the frequency of $x_j$ as predictor of $x_i$ % varFreq[Predictors[i]]+=1
        \end{enumerate}
 \end{enumerate} %end of "repeat the following process"
\item For all other variable $x_j, j=$\{$1,\cdots,p$\}$\setminus$\{$i$\}, if the frequency $x_j$ as predictor of $x_i$ is greater than threshold, set $x_j$ as predictor of $x_i$ in filtered model (one dependency network) %varFreq[]
\end{enumerate}
%=======
%in gibbsreg.cpp
%void gibbsSamplingReg(int targetGene,int* varFreq)
%if(1==VarInModel[Predictors[i]])
%	    {
%               varFreq[Predictors[i]] += 1;
%	    }
%========

In second step, first generate dependency network which is based the filtered model in step1, the possible predictors for initial model are selected from the final model in step1 (in step1, all other variables [$-i$]=\{$1,\cdots,p$\}$\setminus$\{$i$ \} are possible predictors of $x_i$), in addition, the best model is recorded (in step1, final model is average model); then DAG is generated from the best model. In this step, generate several DAGs by running several times of above procedure. For generating compositional network (DAG) from dependency network, see section 4.1.2 in the paper \cite{Dobra_2004}.\\
%==========
%in startpoint.cpp
%//this graph is for keeping the set of possible predictors 
%LPList predictorGraph = InitGraph(NumberOfGenes);
%ReadGraph(predictorGraph,infile);

%in variousfunctions.cpp
%int gibbsSamplingReg(int targetGene,LPList predictorGraph,int* bestModel)
%bestModel[Predictors[i]] = VarInModel[Predictors[i]];
%==========

In third step, the starting points/models (which are generated in step2) are improved in a procedure (simulated annealing) that converges to local models of the posterior distribution over the space of DAGs. \\

After all these three steps, an average network is generated from all statistically significant models during the simulated annealing procedure. Find the model with best posterior probability $w$, all other recorded (for example, recorde every 1000 iterations) models with porstiror probability $w_i$
\begin{equation}
\exp{(w_i-w)}>\epsilon
\end{equation}
will be considered in averaging network.

%\addcontentsline{toc}{chapter}{\numberline{}Bibliography}
%\bibliography{GGM}
\begin{thebibliography}{}

  \bibitem{Dobra_2004}
     Adrian Dobran, Chris Hans, Beatrix Jones, Joseph R. Nevins, Guang Yao \& Mike Wes. 
     {Sparse graphical models for exploring gene expression data}
     \textit{Journal of Multivariate Analysis} \textbf{90} (2004), 196--212.
\end{thebibliography}
\end{document}

