I have multiple Nx3 points, and I sequentially generate a new value for each from its corresponding multivariate Gaussian, each with 1×3 mean and 3×3 cov. So, together, I have arrays: Nx3 array of points, Nx3 array of means and Nx3x3 array of covs.

I only see how to do it with the classic for-loop:

import numpy as np
from scipy.stats import multivariate_normal

# Generate example data
N = 5  # Small number for minimal example, can be increased for real use case
points = np.random.rand(N, 3)
means = np.random.rand(N, 3)
covs = np.array([np.eye(3) for _ in range(N)])  # Identity matrices as example covariances

# Initialize an array to store the PDF values
pdf_values = np.zeros(N)

# Loop over each point, mean, and covariance matrix
for i in range(N):
    pdf_values[i] = multivariate_normal.pdf(points[i], mean=means[i], cov=covs[i])

print("Points:n", points)
print("Means:n", means)
print("Covariances:n", covs)
print("PDF Values:n", pdf_values)

Is there any way to speed this up? I tried to pass all directly to multivariate_normal.pdf, but also from the docs that does not seem to be supported (unlike a simpler case of generating values for Nx3 points, but with same mean and covariance.

Maybe some implementation not from scipy?

I might be too hopeful, but somehow I hope there is a simpler way to speed this up and avoid iterating with this for-loop in a big array of data directly with Pythonic loop.