Multivariate model not behaving as expected

Continuing the discussion from Multivariate model not behaving as expected:

Thanks for your answer! ‘hads’ records the periods and absolute magnitudes of each star (this is a astrophysical project). Actually, if I write the code like this,

likelihood = {}
counter = 0
for i in range(len(band_list)):
    for j in range(len(hads)):
        likelihood[f"m_{counter}"] = pm.Normal(f'm_{counter}',
            mu = mus[j] + alpha[i] * hads['fundamentalized_log_periods'][j] + M_0[i] 
                    + EBV[j] * (extinction_ccm_a[i] * const_R_V + extinction_ccm_b[i]),
            sigma = sigma,
            observed = observed_apparent_mags[counter])
        counter += 1

the result will be alright (but when I use Macbook, it would report ‘bracket nesting level exceeded maximum of 256’), so I think there’s something wrong with the pm.dot. Though the code work in Windows in a dataset of 15 stars, if I expand the dataset to over 100 stars, it would exceed the maximun level. So I may have to use matrix dot product. How could I change the code?