clear all;
% This is an example of using a quasi-Bayesian neighborhood selection
% method for fitting Gaussian graphical models with an elastic net prior.
% For more detail on the methods see the paper
% "A scalable quasi-Bayesian framework for Gaussian graphical models", arXiv:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This code requires Matlab Statistics and Parallel toolbox.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isempty(getenv('PBS_NP'))
NP = 45;
profilename = 'local';
else
NP = str2double(getenv('PBS_NP'));
profilename = 'current';
end
myPool = parpool(profilename, NP);
if isempty(myPool)
error('pctexample:backslashbench:poolClosed', ...
['This example requires a parallel pool. ' ...
'Manually start a pool using the parpool command or set ' ...
'your parallel preferences to automatically start a pool.']);
exit
end
% load cancer group data
% the R code Data_processing.R contains how dtcancer.mat and dtnocancer.mat
% are generated
load 'dtcancer.mat';
sigmasq = []; %we do not know the diagonal. It is then estimated by empirical Bayes
gamma_0 = 0.2;
Niter = 10000;
[sigmasq_est_cancer, CI_L_cancer, CI_m_cancer, CI_R_cancer, Conv_diag_cancer]...
= Bayes_Neigh_Sel(dtcancer, gamma_0, Niter, sigmasq);
% load no-cancer group data
load 'dtnocancer.mat';
[sigmasq_est_no_cancer, CI_L_no_cancer, CI_m_no_cancer, CI_R_no_cancer, Conv_diag_no_cancer]...
= Bayes_Neigh_Sel(dtnocancer, gamma_0, Niter, sigmasq);
%Output:
% sigmasq_est: p by 1: inverse of the diagonal elements of the precision
% matrix
% CI_m: p by p, estimated (sparse) precision matrix (with zero diagonal).
% The diagonal is estimated separately in sigmasq_est
% CI_L and CI_R: p by p, boundaries of the 95% posterior intervals.
% Conv_diag: an integer. The number of neigh selection regression MCMC
% where Geweke convergence fail.
save res_real_data;
%%For post-processing of the results see Data_processing.m
delete(myPool);