Model multivariate normal with separate means, dimension mismatch error

[aseyboldt]

This should do the trick:

mu = tt.stack([mu1, mu2]).T
joint_obs = pm.MvNormal('joint', mu=mu, cov=cov, observed=y)

You might also want to use pm.LKJCholeskyCov instead of pm.LKJCorr, it is usually way faster.