function [lmarglik,b1,s12]=bayesreg(y,x,b0,capv0inv,v0,s02,capv0)

%For a given data set and values for prior hyperparameters
%Calculate OLS and posterior quantities, marginal likelihood
%and predictive (Predictive sets x-star=.5

%Ordinary least squares quantities
k=cols(x);
n=rows(x);
bols = invpd(x'*x)*x'*y;
s2 = (y-x*bols)'*(y-x*bols)/(n-k);
bolscov = s2*invpd(x'*x);
bolssd=zeros(k,1);
for i = 1:k
bolssd(i,1)=sqrt(bolscov(i,i));
end
v=n-k;

%Posterior hyperparameters for Normal-Gamma
xsquare=x'*x;
v1=v0+n;
capv1inv = capv0inv+ xsquare;
capv1=invpd(capv1inv);
b1 = capv1*(capv0inv*b0 + xsquare*bols);
if det(capv0inv)>0
    v1s12 = v0*s02 + v*s2 + (bols-b0)'*invpd(capv0 + inv(xsquare))*(bols-b0);
else
    v1s12 = v0*s02 + v*s2;
end
s12 = v1s12/v1;

bcov = capv1*v1s12/(v1-2);
bsd=zeros(k,1);
for i = 1:k
bsd(i,1)=sqrt(bcov(i,i));
end




%posterior mean and variance of error precision
hmean = 1/s12;
hvar=2/(v1s12);
hsd=sqrt(hvar);


%log of marginal likelihood for the model if prior is informative
if det(capv0inv)>0;
    intcon=gammaln(.5*v1) + .5*v0*log(v0*s02)- gammaln(.5*v0) -.5*n*log(pi);
    lmarglik=intcon + .5*log(det(capv1)/det(capv0)) - .5*v1*log(v1s12);
end