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./chi2rnd(n-k,num,1));



% simulate regression vector



a=chol(RI);

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









