[beta,s]=bay_reg(y,X,num)

function [beta,s]=bay_reg(y,X,num)
%
% BAY_REG Simulated sample from the posterior for a normal regression model.
%	[BETA,S]=BAY_REG(Y,X,NUM) returns a matrix BETA of simulated values
%	from the posterior distribution of the regression parameter and a
%	vector S of simulated values from the posterior of the residual
%	standard deviation, where Y is a column of the observed responses,
%	X is the design matrix, and NUM is the size of the simulated sample.

if nargin==2, num=1000; end
[n,k]=size(X);

% find least-squares estimates, variance matrix, and estimate of residual variance

RI=inv(X'*X);
b = X\y;
s2=1/(n-k)*(y-X*b)'*(y-X*b);

% simulate residual standard deviation

s=sqrt((n-k)*s2./ch2rnd(n-k,num));

% simulate regression vector

a=chol(RI);
beta=randn(num,k)*a.*(s*ones(1,k))+ones(num,1)*b';



function rn=ch2rnd(alpha,n)
%  rn=chi2rnd(alpha,n)generates a vector
%  of n chi2(alpha) variates

if nargin==1, n=1; end

rn=2*rgam(n,alpha/2);

function rn=rgam(n,alpha)
%  rn=rgam(n,alpha)generates a vector
%  of n gamma(alpha) variates
a=alpha-1;
rn=zeros(n,1);
while prod(rn)==0
   v1=-log(rand(n,1));
   v2=-log(rand(n,1));
	  id=v2>=(a*(v1-log(v1)-1));
			rn=rn+v1.*id.*(rn==0);	
end
rn=rn*alpha;