Hi, we are trying to estimate a multivariate normal model: X\sim MVN(\mu,\Sigma). It is understandable that when \mu is passed as an n-dimentional vector, then \Sigma must be n\times n and its element ij denotes covariance between X_i and X_j.

From docs:

cov = np.array([[1., 0.5], [0.5, 2]])

mu = np.zeros(2)

vals = pm.MvNormal(‘vals’, mu=mu, cov=cov, shape=(5, 2))

However, we discovered by chance that \mu does not have to be a vector, but also a squared matrix. The model works, no error message is displayed, and we get quite nice results. But we do not understand HOW the model works exactly.

Questions:

- Does it interpret each row/column of the matrix as a multivariate normal variable?
- And what is the interpetation of \Sigma in this case?

Thanks in advance for any clues

Struggling researchers - Agata and Tomek