I’m a little confused. In your snippet, it appears that chol will only be one 2x2 covariance matrix. I’m not exactly sure what MvNormal.dist is doing; more specifically, I’m not sure what the ‘.dist’ component does? And why pm.Potential?
In the GLM hierarchical linear regression example in the documentation, alpha and beta are both distributed as normal dists with shape = n_countries, which returns different alpha and beta values for each country. Why is something similar not done with pm.LKJCholeskyCov as written in your example? I see you’ve added shape=(n_obs, 2) to pm.MvNormal.dist. What does this function return exactly? When I try running this code, it tells me mvn_dist.shape is equal to [n_obs, 2].