У меня есть несколько точек Nx3, и я последовательно генерирую новое значение для каждой из соответствующего многомерного гауссиана, каждая со средним значением 1x3 и cov 3x3. Итак, у меня есть массивы: массив точек Nx3, массив средних Nx3 и массив covs Nx3x3.
Я вижу только, как это сделать с помощью классического цикла for:
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)
Есть ли способ ускорить это? Я пытался передать все непосредственно в multivariate_normal.pdf, а также из документации, которая, похоже, не поддерживается (в отличие от более простого случая генерации значений для точек Nx3, но с тем же средним значением и ковариацией.
Может быть какая-то реализация не от scipy?
Возможно, я слишком надеялся, но почему-то я надеюсь, что есть более простой способ ускорить это и избежать итерации с этим циклом for в большом массиве данных непосредственно с помощью цикла Pythonic.
🤔 А знаете ли вы, что...
Python поддерживает динамическую типизацию, что облегчает разработку.
Вот решение, включающее разложение Холецкого.
метод 1:
import numpy as np
x = points - means
p = x.shape[1]
res = np.exp(-0.5*(x*np.linalg.solve(covs, x)).sum(1))
res/(2*np.pi)**(p/2)/np.linalg.det(covs)**0.5
array([0.03356053, 0.03167125, 0.08042358, 0.04351325, 0.1328082 ])
метод 2:
y = np.linalg.solve(LU := np.linalg.cholesky(covs), points - means)
np.exp(-0.5* (y**2+np.log(2*np.pi)).sum(1))/np.einsum('kii->ki',LU).prod(1)
array([0.03356053, 0.03167125, 0.08042358, 0.04351325, 0.1328082 ])
метод 3:
multivariate_normal.pdf(y, np.zeros(p))/np.einsum('kii->ki',LU).prod(1)
array([0.03356053, 0.03167125, 0.08042358, 0.04351325, 0.1328082 ])
метод 4: Использование svd
U,S,V = np.linalg.svd(covs)
res_svd = np.einsum('ki,kij,kj,kjl,kl->k', x, U, 1/S, V, x)
np.exp(-0.5 * (res_svd + np.log(2 * np.pi * S).sum(1)))
array([0.03356053, 0.03167125, 0.08042358, 0.04351325, 0.1328082 ])
Сгенерировать данные:
import numpy as np
from scipy.stats import multivariate_normal, wishart
# Generate example data
N = 5 # Small number for minimal example, can be increased for real use case
p = 3 # dimension
np.random.seed(10)
points = np.random.rand(N, p)
means = np.random.rand(N, p)
covs = wishart.rvs(p, np.eye(p),N, 2) # 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])
pdf_values
array([0.03356053, 0.03167125, 0.08042358, 0.04351325, 0.1328082 ])