About set of Multivariate Normal Distributions

Questions
1

Let’s say we have a matrix:

M = \begin{bmatrix} \mu_{11} & \mu_{12} & \dots & \mu_{1m} \\ \mu_{21} & \mu_{22} & \dots & \mu_{2m} \\ \hdotsfor{4} \\ \mu_{n1} & \mu_{n2} & \dots & \mu_{nm} \end{bmatrix}

and a matrix:

\Sigma = \begin{bmatrix} \sigma_{11} & \sigma_{12} & \dots & \sigma_{1m} \\ \sigma_{21} & \sigma_{22} & \dots & \sigma_{2m} \\ \hdotsfor{4} \\ \sigma_{m1} & \sigma_{m2} & \dots & \sigma_{mm} \end{bmatrix}

where each row \inline (\mu_i_1, \mu_i_2, \dots, \mu_i_m) of a matrix \inline M is supposed to be a mean vector of some m-dimensional Multivariate Normal Distribution with Covariance matrix \inline \Sigma.

And we need to define \inilin n Multivariate Normal Distributions with different mean vectors but with the same covariance matrix inside of a Model. How to achieve this without explicit "for loop" which is a quite performance overhead?

A rough code:

with pm.Model():
    Sigma = pm.Deterministic(...)
    M = pm.Deterministic(...)
    
    # unoptimal solution
    for i in xrange(n):
        _Y = pm.MvNormal('_Y_' + str(i), mu=mu_[:, i], cov=sd_, shape=(m, ))
    
    # Is it possible to vectorize last lines using pm.MvNormal or pm.MatrixNormal ??

Thanks,
Best Regards
Alex

1 Like
2

If the MvNormal is a latent variable, it is much more straightforward to reparameterized it as Y = mu + L.dot(s) with s ~ Normal(0, 1)

n, m = 100, 5
Mean = np.random.randn(m, n)*5

with pm.Model() as model:
    sd_dist = pm.HalfNormal.dist(2.5)
    packed_chol = pm.LKJCholeskyCov('chol', n=m, eta=2, sd_dist=sd_dist)
    # compute the covariance matrix
    L = pm.expand_packed_triangular(m, packed_chol, lower=True)
    mu = pm.Normal('mu', 0., 1., shape=(m, n))
    Mv = pm.Deterministic('Mv', L.dot(mu) + Mean)

Using MvNormal directly is fine too:

n, m = 100, 5
Mean = np.random.randn(m, n)*5
with pm.Model() as model:
    sd_dist = pm.HalfNormal.dist(2.5)
    packed_chol = pm.LKJCholeskyCov('chol', n=m, eta=2, sd_dist=sd_dist)
    # compute the covariance matrix
    L = pm.expand_packed_triangular(m, packed_chol, lower=True)
    Mv = pm.MvNormal('Mv', Mean.T, chol=L, shape=(n, m))

Note that the shape from MvNormal is different compare to above where you use dot product.

3 Likes
3

Thanks! The second solution is exactly what I need.

Best Regards
Alex