Choose out of multivariate normal pdf

I have some questions about scipy.stats.multivariate_normal.

1. What does it return? I understand that when using:

   var = scipy.stats.multivariate_normal(mean, cov)
   z = var.pdf([x,y])
   

   I'll get in z the pdf of this index. Is this true?

2. I assumed that the mean given is the peak of the pdf. Am I right?

3. I want to create a 4D map of gaussian pdf, and choose out of it randomly, by the pdf, N points (for particle filter predicting). To use np.choose I need to create an array of the pdf values of shape (M,1), and an array of the 4D points of shape (M, 4).

Is there any more efficient way to do this?

The normal distribution is quite ubiquitous and implemented in many different places, so I think you would be better off using the random normal sampling from numpy: https://numpy.org/doc/stable/reference/random/generated/numpy.random.normal.html. In this function you can specify the shape of the output.