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