%clear all;

% This is an example of using the Moreau posterior approx. 
% for empirical Bayes inference of linear regression models with an elastic net prior. 
% For more detail on the method see the paper
% "A Moreau-Yosida approximation scheme for posterior quasi-posterior distributions", arXiv:1505.07072. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This code requires Matlab Statistics toolbox.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% Linear regression example
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Data input should be:  y = n_by_1 vector, and X = n_by_p

y = dtcancer(:,10);
X = dtcancer(:,[1:9,11:984]);
Niter = 1000; %MCMC iterations
gamma_0 = 0.2; % best choice between 0.1 and 0.25
sigmasq_star = []; %if [] then empirical Bayes is used

[est_delta, Res_CI, Res_others, Res_diag] ...
                       = Moreau_Approx_lm(y, X, gamma_0, Niter, sigmasq_star);


%%Output:
%est_delta: estimated model structure and probability of equal to 1, a p by 2 matrix.
%Rest_CI: 95% posterior confidence intervals, and the posterior mean. A p by 3 matrix.
% Res_others: MCMC output for lambda_1, lambda_2 and q (see the paper for explanation).
% Res_diag: Diagnostic output, a 5 by 1 vector. Contains:
%    (*) Geweke convergence diagnostic output (1= 'we do not reject stationarity' at 10%), 
%    (*) Acceptance probab. of the three type of MCMC moves3 average,
%    (*) and the estimated model noise variance sigma_sq.



