gofg_sp.m

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

function [r, Qqr, Rqr]= ...
gofg_sp(gamma,indice,X,Ytilde,QqrP,RqrP,elem,k,delta,w,v1);

%gofg_sp: Selection prior - Function that
%         computes the log-relative posterior probability
%*******************************************************************
%INPUTS:
% gamma,          binary vector
% indice,         1 -- QR decomp
%                 2 -- QRinsert
%                 3 -- QRdelete
% X, 		  independent variables, n by p
% Ytilde,         augmented response matrix, (n+p) by q
% QqrP,RqrP,      QR-matrices of the previous vector
% elem,           position to add/delete
% k,delta,        hyperparameters Inverse Wishart
% w,              hyperparameter Bernoulli priors
% v1,             hyperparameter - normal selection prior
%                 (column vector (p by 1) of standard deviations)
%
%OUTPUTS:
% r,              log-relative posterior probability of gamma
% Qqr, Rqr,       QR matrices of gamma
%
%USAGE:
% [r, Qqr, Rqr]= ...
% gofg_sp(gamma,indice,X,Ytilde,QqrP,RqrP,elem,k,delta,w,v1);
%
%NOTES:
% Called by bvsgs_sp and itergs_sp
%
%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);
Qk=k*eye(q);
I=find(gamma);
nonnulle=length(I);
switch indice
        case 1  % ------------------ Indice = 1 -- QR Decomposition
		sotto=zeros(p,nonnulle);
		sotto(I,:)=eye(nonnulle);
			% Augmented X, (nOss+p) by nonnulle, 
			% H0=Diag*R*Diag,R=Id
		Xtilde=cat(1,X(:,I)*diag(v1(I)), sotto); 
 		[Qqr Rqr]=qr(Xtilde);        

        case 2  % ------------------ Indice = 2 -- QRinsert
		prima=find(I<elem);
		jj=length(prima)+1;
		sotto=zeros(p,1);
		sotto(elem)=1;
		ins=cat(1, X(:,elem)*v1(elem), sotto);
		[Qqr Rqr]=qrinsert(QqrP,RqrP,jj,ins);

                                
        case 3  % ------------------ Indice = 3 -- QRdelete
		prima=find(I<elem);
		jj=length(prima)+1;
		[Qqr Rqr]=qrdelete(QqrP,RqrP,jj);      
end
sub=Qqr(:,nonnulle+1:nOss+p)'*Ytilde;	% Seber, page 498
RSP=sub'*sub;                       	% residual sum of products
Qstar=Qk+RSP;                       	% q by q
switch nonnulle
	case 1
        	detH0K=Rqr(1,1)^2;
	otherwise
        	detH0K=prod(diag(Rqr).^2);
end
prog=sum(log((w.^gamma).*(1-w).^(1-gamma)));
r= -(q/2)*log(detH0K) - ((delta+nOss+q-1)/2)*log(det(Qstar)) + prog;