function [sigmasq_est, CI_L, CI_m, CI_R, Conv_diag] = Bayes_Neigh_Sel(X, gamma_0, Niter, sigmasq)
%Moreau_Approx_ggm(X, gamma_0, Niter, sigmasq, theta_star)

p = length(X(1,:));

CI_L = zeros(p-1,p);
CI_m = zeros(p-1,p);
CI_R = zeros(p-1,p);
theta=zeros(p-1,p); 
struc_theta=zeros(p-1,p);

Geweke = zeros(1,p);

sigmasq_est=zeros(p,1);

parfor jj = 1:p
    if jj == 1
        dep_var = 2:1:p;
    elseif jj==p
        dep_var = 1:1:(p-1);
    else
        dep_var = [1:1:(jj-1) (jj+1):1:p];
    end
    if isempty(sigmasq)
        sigmasq_j=[];
    else
        sigmasq_j=sigmasq(jj);
    end
    [dt, Res_CI, ~, Res_diagnostic]  ...
             = Moreau_Approx_lm(X(:,jj), X(:,dep_var), gamma_0, Niter, sigmasq_j);
         
    sigmasq_est(jj) = Res_diagnostic(5);
    theta(:,jj) = -(1/sigmasq_est(jj))*Res_CI(:,2);
    struc_theta(:,jj) = dt(:,1);
    CI_L(:,jj) = -(1/sigmasq_est(jj))*Res_CI(:,3); 
    CI_m(:,jj) = -(1/sigmasq_est(jj))*Res_CI(:,2); 
    CI_R(:,jj) = -(1/sigmasq_est(jj))*Res_CI(:,1);
    Geweke(jj) = 1 - Res_diagnostic(1);
end

Conv_diag = sum(Geweke);

%estimate parameter and graph structure
theta = [zeros(1,p);theta];
struc_theta = [zeros(1,p);struc_theta];
CI_L = [zeros(1,p);CI_L];
CI_m = [zeros(1,p);CI_m];
CI_R = [zeros(1,p);CI_R];


for j =1:p
    if j == 1
        dep_var = 2:1:p;
    elseif j==p
        dep_var = 1:1:(p-1);
    else
        dep_var = [1:1:(j-1) (j+1):1:p];
    end
    theta(dep_var,j) = theta(2:p,j);
    theta(j,j) = 1/sigmasq_est(j);
    struc_theta(dep_var,j) = struc_theta(2:p,j);
    struc_theta(j,j)=1;
    CI_R(dep_var,j)= CI_R(2:p,j);
    CI_m(dep_var,j)= CI_m(2:p,j);
    CI_L(dep_var,j)= CI_L(2:p,j);
end

struc_theta(find(struc_theta>=0.5))= 1;
struc_theta(find(struc_theta<0.5))= 0;
struc_theta = struc_theta.*(struc_theta'); %symmetrization
theta = 0.5*(theta + theta');
theta = theta.*struc_theta;


%Confidence  intervals
CI_L = bsxfun(@min,CI_L,CI_L'); CI_m = 0.5*bsxfun(@plus,CI_m,CI_m');
CI_R = bsxfun(@max,CI_R,CI_R');
CI_L = CI_L.*struc_theta;
CI_m = CI_m.*struc_theta;
CI_R = CI_R.*struc_theta;

end
