As another update, replacing:
betas = at.reshape(at.tile(alphas.T,num_respondents), (num_respondents, num_items)) + at.dot(betas_init,L_sigma)
with:
betas = alphas.reshape(num_respondents, 1) + at.dot(betas_init,L_sigma)
Allowed me to sample from a standard normal.