function [tau2,ialpha_prvar] = shrink_IG(alphadraw,p,plag,phi1,phi2)
   
a1 = 0.5 + phi1;
a2 = 0.5*(alphadraw.^2) + phi2*[ones(p,1) ; kron((1./(1:plag).^2)',ones(p*p,1))];
tau2 = 1./gamrnd(a1,1./a2);

% Restrict tau2 for numerical stability
tau2(tau2>=2) = 2;
tau2(tau2<=1e-10) = 1e-10;

D0 = diag(tau2);
D0(1:p,1:p) = 2*eye(p);  % Intercepts are unrestricted

ialpha_prvar = inv(D0);