Python: Sample from multivariate normal with N means and same covariance matrix

Question

Suppose I want to sample 10 times from multiple normal distributions with the same covariance matrix (identity) but different means, which are stored as rows of the following matrix:

means = np.array([[1, 5, 2],
                  [6, 2, 7],
                  [1, 8, 2]])

How can I do that in the most efficient way possible (i.e. avoiding loops)

I tried like this:

scipy.stats.multivariate_normal(means, np.eye(2)).rvs(10)

and

np.random.multivariate_normal(means, np.eye(2))

But they throw an error saying mean should be 1D.

Slow Example

import scipy
np.r_[[scipy.stats.multivariate_normal(means[i, :], np.eye(3)).rvs() for i in range(len(means))]]

in this special case, isn't one sample is like adding this means matrix to a 3x3 normal N(0,1) matrix? — Lior Cohen, Dec 11 '20 at 14:08
How about we start with - what shape of output do you expect? Because it's not at all clear what you want to get out of this. — Daniel F, Dec 11 '20 at 14:13

Quang Hoang · Answer 1 · 2020-12-11T15:44:44.540

1

Your covariance matrix indicate that the sample are independent. You can just sample them at once:

num_samples = 10
flat_means = means.ravel()

# build block covariance matrix
cov = np.eye(3)
block_cov = np.kron(np.eye(3), cov)

out = np.random.multivariate_normal(flat_means, cov=block_cov, size=num_samples)

out = out.reshape((-1,) + means.shape)

edited Dec 11 '20 at 15:44

answered Dec 11 '20 at 15:34

Quang Hoang

146,074
10
56
74

They are independent only cause this is a MWE. I want it to work with a full covariance matrix. Sorry for the confusion – Euler_Salter Dec 11 '20 at 15:35
@Euler_Salter use block matrix to construct the `9x9` covariance matrix. – Quang Hoang Dec 11 '20 at 15:36
But they won't be independent – Euler_Salter Dec 11 '20 at 15:37
Oh I see what you mean. Surely though the introduction of all those extra entries will slow it down compared to a list comprehension? – Euler_Salter Dec 11 '20 at 15:43
It might when you draw small samples. But I bet it would be much faster when you need to draw more samples. Try with 1000 samples and you see that this would be hundred times faster. – Quang Hoang Dec 11 '20 at 15:45

Python: Sample from multivariate normal with N means and same covariance matrix

Slow Example

1 Answers1

Linked

Related