bvsgs_sp.m

Bayes model averaging with selection of regressors - This program can be utilized for Bayesian Varia

function [Gamma, GammaD, logProb, logProbD, PostGamD, MargGam, SWITCH]= ...
bvsgs_sp(gamprec,X,Y, delta,k, w, v1)

%bvsgs_sp: Bayesian Variable Selection - Gibbs Sampler - Main Program
%	   Selection prior (that is, v_0j = 0, j=1...p)
%	   Bernoulli priors with w_j = w, j=1...p
%***********************************************************************
%INPUTS:
% gamprec,      starting binary vector.
%               If gamprec=[] the program asks for r to select a starting
%               vector with first r elements equal to 1 
% X,            independent variables, matrix (n by p)
% Y,            response variables, matrix (n by q)
% delta, k,     hyperparameters Inverse Wishart
% w,   		Bernoulli prior 
% v1,		hyperparameter - normal selection prior 
%		(column vector (p by 1) of standard deviations)
%
%OUTPUTS:
% Gamma,        all visited vectors - matrix in sparse form
%		(Gamma(1,:) contains the starting vector)
% GammaD,	distinct visited vectors, ordered according to
%		their (normalized) relative posterior probability
%		(matrix in sparse form)
% logProb,      log-relative posterior probabilities 
%		of all visited vectors
% logProbD	log-relative posterior probabilities
%               of distinct visited vectors
% PostGamD	normalized ordered relative probabilities 
%		of distinct visited vectors 
% MargGam	marginal probabilities of components 
% SWITCH,	number of component switches (out of p) from iteration 
%		to iteration.
%		
%USAGE:
% [Gamma, GammD, logProb, logProbD, PostGamD, MargGam, SWITCH]= ...
%	bvsgs_sp(gamprec, X,Y, delta,k, w, v1)
%
%NOTES:
% Data must be centered.
% The programs asks for possible permuting of the data and for the
% Gibbs parameters (initial number of variables included, number
% of iterations).
% QR matrices updated every m iterations (m provided by the user).
% Gamma and GammaD are in sparse form.
% Functions called: gofg_sp, gibbs_sp, repliche, probord
%
%REFERENCE: 
% Brown, P.J., Vannucci, M. and Fearn, T.
% Multivariate Bayesian variable selection and prediction
% Journal of the Royal Statistical Society B, 60(3), 1998, pp. 627-641. 
%
%Copyright (c) 1997 Marina Vannucci.
%**********************************************************************

[nOss p]=size(X);
[q]=size(Y,2); 
                % p independent variables, q responses, nOss observations
disp(' ')
disp('------- GIBBS SAMPLER --------')
disp(' ')

domanda= ...
input('-- Do you want to permute the order of the X variables?(y/n)','s');
switch domanda
	case 'y'
		indici=randperm(p);
		X=X(:,indici);
	case 'n'
		indici=1:p;
end

if isempty(gamprec)
	disp(' ')
	disp('-- Starting binary vector with first r elements equal to 1');
	r=input('   please provide r, r=?');	
	gamprec=[ones(1,r) zeros(1,p-r)];
else
	gamprec=gamprec(indici);
end

disp(' ')
disp('-- To avoid build up of rounding errors,') 
disp('   QR matrices are recomputed every m iterations')
nPar=input('   please provide  m, m=?'); 
disp(' ')
ngibbs=input('-- Number ot total iterations (multiple than m)=?');

disp(' ')
disp('------- RUNNING GIBBS SAMPLER ...')
disp(' ')

v1=v1(indici);       	% standard deviations (mixture normal prior)
Gamma=sparse(zeros(ngibbs+1, p));
logProb=zeros(1,ngibbs+1);
Gamma(1,indici)=gamprec;
SWITCH=zeros(1,ngibbs);
conta=1;
Ytilde=cat(1,Y, zeros(p,q));      % augmented response matrix

while conta<ngibbs+1

%-------------------- updating QR decomposition
	[probprec, Qqr, Rqr] = ...
		gofg_sp(gamprec,1,X,Ytilde,[],[],[],k,delta,w,v1);
	logProb(conta)= probprec;

%-------------------- Gibbs (nPar iterations)
        [ProbSt, GammaSt, storeS]= gibbs_sp( ...
	gamprec,probprec,X,Ytilde,Qqr,Rqr,k,delta,w,v1,nPar);
	Gamma(conta+1:conta+nPar,indici)=GammaSt;
	logProb(conta+1:conta+nPar)=ProbSt;
	SWITCH(conta:conta+nPar-1)=storeS;
	conta=conta+nPar;
	gamprec=GammaSt(nPar,:);
	disp(' ')
        disp('Number of iterations:')
        disp(conta-1)
        disp('...')
end

disp(' ')
disp('------- Searching for replicates ...')
disp(' ')
[GammaD, logProbD]=repliche(Gamma, logProb);

%------------------- Normalized posterior and marginal probs.
%		     Ordering of distinct visited vectors according
%                    to probability

[GammaD, logProbD, PostGamD, MargGam]=probord(logProbD, GammaD);

disp(' ')
disp('--- End.')