function h = drawh(y,X,rhoVec,m,v)
%{
draw precision
modelled after Jeremy Piger's gen_h.m

Dick Startz
May 2016
modified September 2016
modified March 2017 to allow vector of rho's
%}

e = y - X*rhoVec;
vbar = size(e,1) + v;
mbar = vbar/(e'*e + v/m);

%{
 The following converts mbar and vbar into the definition of the gamma
 definition used by Matlab random number generator (gamrnd), Gamma(a,b).
 In this paramterization of Gamma, E(h) = a*b and Var(h) = a*b^2.
 a is called the shape parameter and b the scale parameter.
%}
a = vbar/2;
b = mbar/a;

h = gamrnd(a,b); %h is a draw from the posterior for the error precision
end

