function [sigma,naccept2]=draw_sigma(M,iscale,T1,resid,Qstar,beta2m,sigmam,sigmaold)
naccept2=0;
sigma=sigmaold;
sigmanew=iwpq(T1,iscale);
A0new=chol(sigmanew,'lower');
numerator=loglik_m_given_y(M, resid, A0new, Qstar', beta2m, sigmam.^0.5);

A0old=chol(sigmaold,'lower');
denominator=loglik_m_given_y(M, resid, A0old, Qstar', beta2m, sigmam.^0.5);

accept=min([exp(numerator-denominator) 1]);
u=rand(1,1);
if u<accept
sigma=sigmanew;
naccept2=1;
end

