function [beta,v,prob,haty,e,mliktau]=mregfullbayes(y,x,tau,A)
% fits a multiple linear regression model to response data y 
%  using predictor matrix x 
%
% Conjugate shrinkage prior: zero mean, common variance tau for all predictors
% 	note that the intercept has a very large prior variance -- vague on intercept
% 	Residual precision has Gamma(a/2, b/2) prior : A=[a b]
% Also computes model marhinal likelihood -- for range of possible tau values in Tau
%
%  ouputs are the LSE beta, the SDs of elements of beta in V, 
%  approximate probability level for each coeff, fitted values haty and residuals e
%
%  y is nx1 column vector, x is pxn matrix 
%  plots betas and approximate 90% intervals


[p,n]=size(x); X=x; 		

% model fit ...
H=[ones(n,1),X'];       	        % create design matrix with intercept and center genes
C=eye(p+1)/tau; C(1,1)=1e-12;	 	 % prior precision mx
V=inv(H'*H+C); beta=V*(H'*y);	  	 % compute posterior mean/var
					 % n.b. note its a trivial change to LSE 

haty=H*beta; 				% fitted values
e=y-haty;  				% residuals

ag=A(1); bg=A(2); 	 		% prior for residual variance - deals with proper and ref
if (ag==0 | bg==0) 
	ag=0; bg=0; k=0;
else
	k=(ag/2)*log(bg/2)-gammaln(ag/2); 	 
end
	
ag=ag+n; bg=bg+y'*e; 			% posterior for residual variance
s=bg/ag;		 		% posterior mean for residual precision

mliktau=k+gammaln(ag/2)+0.5*(log(det(C))+log(det(V))-ag*log(bg));
					% log of model marginal likelihood 

v=sqrt(s*diag(V));              	% compute standard errors of beta coefficients
df=n;
prob=2*(1-pt(abs(beta./v),df));         % prob level for each coeff 
qtc=qt(0.95,df); 			% 90% intervals





figure(1) 
plot([0 1+p],[0 0]); hold on
for j=1:p
   plot([j,j],[beta(j+1)+qtc*v(j+1),beta(j+1)-qtc*v(j+1)],'k')
   text(j,beta(j+1),'+','HorizontalAlignment','Center')
end;
xlabel('Predictor gene index'); ylabel('Coeff')
title('Regression coefficients') 
hold off;