Well if you have a station where you can write the looped and non looped version together in a model (probably a smaller version of the data you are using to make things easier), it is easy to debug where the two values deviate from each other using .eval() statement on random variables. Just be careful that if you say for instance RV.eval() and RV depends on some other parameter a, if you do not initialize a then it will randomly assign it so better to say something like RV.eval({a:1}) (note a is supplied to the dict as a variable not as a string name).
However if i’s are your rows and j’s columns then your normal’s mu should be roughly be
mus[None,:] + alpha[:,None]*Lp[None,:] +
M_0[:,None] +
EBV[None,:]*(ExtA[:,None] * RV + ExtB[:,None])
where Lp = logperiods, ExtA = extinction_ccm_a etc. And then you can use observed as it is modulo correctly specifying the shape of the likelihood maybe. Can’t check without a fully running code though.
In general if you want to get two (one dimensional) vectors A, B and want to multiply them to get a matrix whose i,j th element is A_iB_j then you can do M = A[:,None] * B[None,:]. I suspect that is what you need to do here a lot to vectorise things.