Numpy/Scipy: 奇异矩阵计算多变量观测的概率

Question

我正在尝试计算矩阵中观察值的概率，其中我的行是观察值，我的列是使用 python 的特征。我总是遇到奇异矩阵错误，即使在使用随机矩阵时也是如此，所以我怀疑我的代码有问题：

from scipy.stats import multivariate_normal
import numpy as np


def get_mean(x, axis=0):
  return x.mean(axis=axis)


def get_sigma(x):
  return np.cov(x, rowvar=False)


def get_probabilities(x, mu, sigma):
  return multivariate_normal.pdf(x, mean=mu, cov=sigma)


x = np.random.rand(10,10)
t = np.random.rand(1, 10)

mu = get_mean(x)
sigma = get_sigma(x)
p = get_probabilities(t, mu, sigma)

这导致：

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<stdin>", line 2, in get_probabilities
  File "/usr/local/lib/python3.5/dist-packages/scipy/stats/_multivariate.py", line 512, in pdf
    psd = _PSD(cov, allow_singular=allow_singular)
  File "/usr/local/lib/python3.5/dist-packages/scipy/stats/_multivariate.py", line 159, in __init__
    raise np.linalg.LinAlgError('singular matrix')
numpy.linalg.linalg.LinAlgError: singular matrix

我做错了什么？

Answer 1

我认为问题在于 t 是形状 (1,10) 的 np.array 而它应该是形状 (10,)。

如果替换代码行

t = np.random.rand(1, 10)

来自

t = np.random.rand(10)

这有效。

Answer 2

没有足够的数据来估计协方差矩阵 sigma。 x 需要比列更多的行来估计定义明确的协方差矩阵。

例如，这将使代码运行没有问题：

x = np.random.rand(100, 10)

或者，如果您只是没有足够的数据，您可以使用 regularization（或收缩）来更好地调节协方差矩阵：

sigma += np.eye(10) * 1e-3  # problem: how to chose the factor

Scikit-learn 具有 ledoit_wolf 协方差估计器，可根据可用数据自动确定正则化量。

Numpy/Scipy: 奇异矩阵计算多变量观测的概率

Numpy/Scipy: Singular Matrix Calculating probability of multivariate observation

python

numpy

probability

scipy