pbvs_sp.m

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

function [BayesPred,LSPred,ILS,IB]=pbvs_sp(X,Y,Xf,Yf,PostGamD,GammaD,numero,v1)

%pbvs_sp: Prediction - Selection Prior
%*******************************************************************
%INPUTS:
% X,Y,          training data
% Xf,Yf,        validation data
% PostGamD,     normalized ordered relative probabilities
%               of distinct visited vectors
% GammaD,       distinct visited vectors, ordered according to
%               their (normalized) relative post. prob. (PostGamD)
% numero,       number of most likely models for Bayes prediction
% v1            hyperparameter - normal selection prior
%               (column vector (p by 1) of standard deviations)
%
%OUTPUTS:
% BayesPred	Bayes prediction with the 'numero' most likely models
% LSPred        Least Squares Prediction with the best model
% ILS           indices of selected variables for LS prediction
% IB		indices of selected variables for Bayes prediction
%
%USAGE:
% [BayesPred,LSPred,ILS,IB]=pbvs_sp(X,Y,Xf,Yf,PostGamD,GammaD,numero,v1)
%
%NOTES:
% Training Data must be centered
% Validation data must be centered with training means
%
%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);          
[nOssP pP]=size(Yf);
ILS=[];
IB=[];
Ytilde=cat(1,Y, zeros(p,q));
%------------------------------------- Bayes prediction
pred=zeros(nOssP,pP);
costante=sum(PostGamD(1:numero));
for i=1:numero
	II=find(GammaD(i,:));
	nonnulle=length(II);
	H0=diag(v1(II).^2);    
	XH0=X(:,II)*sqrt(H0);
	sotto=zeros(p,nonnulle);
	sotto(II,:)=eye(nonnulle);    
	Xtilde=cat(1, XH0, sotto);
	Xfg=Xf(:,II);
	pred1=Xfg*sqrt(H0)*inv(Xtilde'*Xtilde)*Xtilde'*Ytilde;
	pred=pred+pred1*PostGamD(i)/costante;
	IB=union(IB,II);
end
residual=Yf-pred;
BayesPred=diag(residual'*residual)./nOssP;

%-------------------------------------- LS prediction
ILS=find(GammaD(1,:));
Xfg=Xf(:,ILS);
Xg=X(:,ILS);
residual=Yf-Xfg*inv(Xg'*Xg)*Xg'*Y;
LSPred=diag(residual'*residual)./nOssP;