function [Mbeta,Ms,Mlam]=rreg(y,X,v,m)

%

% RREG  Simulated sample from the posterior for a regression model with t errors.

%	[MBETA,MS,MLAM]=RREG(Y,X,V,M) returns a matrix MBETA of simulated values

%	from the posterior distribution of the regression parameter, a

%	vector MS of simulated values from the posterior of the scale parameter,

%	and a vector MLAM of posterior means of the observation scale parameters, 

%	where Y is a column of the observed responses, X is the design matrix, 

%	V is the degrees of freedom of the t density and M is the size of the simulated 

%	sample.



if nargin<=3, m=1000; end

if nargin<=2, v=10; end

[n,k]=size(X);



Mbeta=zeros(m,k); Ms=zeros(m,1); Mlam=zeros(n,1);



beta=inv(X'*X)*(X'*y);

s2=1/n*(y-X*beta)'*(y-X*beta);

v1=(v+1)/2;



for i=1:m

	lam=gamrnd(v1,2./((y-X*beta).^2/s2+v));

	

	W=diag(lam);

	var=inv(X'*W*X);

	b=var*(X'*W*y);

	a=chol(var);

	beta=sqrt(s2)*a'*randn(k,1)+b;



	s2=1./gamrnd(n/2,2/sum(lam.*(y-X*beta).^2));



	Mbeta(i,:)=beta';

	Ms(i)=sqrt(s2);

	Mlam=Mlam+lam;

end



Mlam=Mlam/m;





