
function [beta, btilde, Btilde] = drawbeta(z,x, w, tau2, theta, invB0, invB0b0)

% Output
%--------------------------------------------------------------------------
% beta      : Gibbs draw of beta, a column vector
% btilde    : posterior mean, a coulmn vector
% Btilde    : posterior covariance matrix
% Input
% z         : Gibbs draw of latent response variable, a column vector
% x         : covariates including a column of ones
% w         : a column vector
% theta     : parameter from normal-exponential mixture of Laplace dist.
% b0        : prior mean
% B0        : prior covariance matrix
%--------------------------------------------------------------------------

[n,k]   = size(x);

 meancomp = zeros(n,k);
 varcomp  = zeros(k,k,n);

for j = 1:n
    meancomp(j,:)  = ( x(j,:).*( z(j) - theta*w(j) ) )./( tau2*w(j) );
    varcomp(:,:,j) = ( x(j,:)'*x(j,:))./( tau2*w(j) );   
end

sum(meancomp);
% posterior variance
Btilde    = ( invB0 + sum(varcomp,3))\eye(k);   % sum(varcomp,3) sums the elements of the kxk matrices

% posterior mean
btilde    = Btilde*( invB0b0 + ( sum(meancomp))' );

% btilde    = Btilde*( B0\b0' + ( sum(meancomp))' );  % Using this slows
% down computation.
% beta      = mvnrnd(btilde',Btilde);

L         = chol(Btilde,'lower');                    % produces a lower triangular matrix
beta      = btilde + (L*mvnrnd(zeros(1,k),eye(k))');   % Gibbs draw

end

