itergs_sp.m

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

function [gamnew, probnew, QqrN, RqrN, swit]= ...
itergs_sp(gamprec,probprec,X,Ytilde,QqrP,RqrP,k,delta,w,v1);

%itergs_sp: Gibbs Sampler - Selection prior - Single iteration
%*******************************************************************
%INPUTS:
% gamprec,        previous binary vector.
% probprec,       log-relative posterior probability of gamprec
% X,              independent variables, n by p
% Ytilde,         augmented response matrix, (n+p) by q
% QqrP,RqrP,      QR-matrices for gamprec
% k,delta,        hyperparameters Inverse Wishart
% w,              hyperparameter Bernoulli prior
% v1,             hyperparameter - normal selection prior
%                 (column vector (p by 1) of standard deviations)
%
%OUTPUTS:
% gamnew,         new visited vector
% probnew,        log-relative posterior probability of gamnew
% QqrN, RqrN,	  QR-matrices for gamnew
% swit,           number of component switches (out of p) 
%		  for this iteration
%
%USAGE:
% [gamnew, probnew, QqrN, RqrN, swit]= ...
% itergs_sp(gamprec,probprec,X,Ytilde,QqrP,RqrP,k,delta,w,v1);
%
%NOTES:
% One iteration consists of p Bernoulli draws.
% Called by gibbs_sp
% Functions called: bernoulli.m, gofg_sp.m
%
%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(Ytilde,2);
                % p independent variables, q responses, nOss observations
swit=0;
%-------------------sampling from the full conditionals of single elements
for j=1:p
	switch gamprec(j)

%---------------case previous j-th element equal 0
		case 0
			probzero=probprec;
			gamtmp=gamprec;
			gamtmp(j)=1;
			[probone, Qqrtmp, Rqrtmp]= ...
			gofg_sp(gamtmp,2,X,Ytilde,QqrP,RqrP,j,k,delta,w,v1); 

%---------------sampling from the Bernoulli distribution
			teta=1/(exp(probzero-probone)+1);
                	gamprec(j)=bernoulli(teta);
                        
%---------------updating posterior probability for the next component 
                	if gamprec(j)==1
                        	probprec=probone;     
                                QqrP=Qqrtmp;
                                RqrP=Rqrtmp;
				swit=swit+1;
			else
                        end
 
%---------------case previous j-th element equal 1
		otherwise
			probone=probprec;
			gamtmp=gamprec;
                        gamtmp(j)=0;
			[probzero, Qqrtmp, Rqrtmp]= ...
			gofg_sp(gamtmp,3,X,Ytilde,QqrP,RqrP,j,k,delta,w,v1);
                
%---------------sampling from the Bernoulli distribution                        
                        teta=1/(exp(probzero-probone)+1);
                        gamprec(j)=bernoulli(teta);

%---------------updating posterior probability for the next component
                        if gamprec(j)==0
                                probprec=probzero;
                                QqrP=Qqrtmp;
                                RqrP=Rqrtmp;
				swit=swit+1;
                        else
			end
	end	
end
gamnew=gamprec;
probnew=probprec;
QqrN=QqrP;
RqrN=RqrP;