我試圖使用 scipy 來估計數據集在某些點的密度。from scipy.stats import gaussian_kdeimport numpy as np我有一個A3D 點數據集（這只是一個最小的例子。我的實際數據有更多維度和更多樣本）A = np.array([[0.078377 , 0.76737392, 0.45038174], [0.65990129, 0.13154658, 0.30770917], [0.46068406, 0.22751313, 0.28122463]])以及我想要估計密度的點B = np.array([[0.40209377, 0.21063273, 0.75885516], [0.91709997, 0.79303252, 0.65156937]])但我似乎無法使用該gaussian_kde功能，因為result = gaussian_kde(A.T)(B.T)回報LinAlgError: Matrix is not positive definite我該如何修復這個錯誤？如何獲得樣品的密度？
1 回答
波斯汪
TA貢獻1811條經驗 獲得超4個贊
您的數據中具有高度相關的特征，這會導致數值錯誤。有幾種可能的方法可以解決這個問題，每種方法都有優點和缺點。下面提出了一個直接替換類gaussian_kde。
診斷
您的（即您在創建對象時而不是dataset使用對象時提供的矩陣）可能包含高度依賴的特征。這一事實（可能與低數值分辨率和“太多”數據點相結合）導致協方差矩陣具有接近零的特征值，由于數值誤差可能會低于零。這反過來會破壞在數據集協方差矩陣上使用 Cholesky 分解的代碼（具體細節請參閱解釋）。gaussian_kdefloat32
假設您的數據集具有形狀，(dims, N)您可以通過 測試這是否是您的問題np.linalg.eigh(np.cov(dataset))[0] <= 0。如果有任何結果出來True，請讓我第一個歡迎您加入俱樂部。
治療方法
這個想法是讓所有特征值都大于零。
將數值分辨率提高到實用的最高浮點數可能是一個簡單的解決方案，值得嘗試，但可能還不夠。
考慮到這是由相關特征引起的，刪除數據點并沒有多大幫助。通過減少需要粉碎的數字來傳播更少的誤差并保持特征值大于零的希望渺茫。它很容易實現，但它會丟棄數據點。
更復雜的修復方法是識別高度相關的特征并將它們合并或忽略“多余”的特征。這是很棘手的，尤其是當維度之間的相關性因實例而異時。
也許最干凈的方法是保持數據不變，并將有問題的特征值提升為正值。這通常通過兩種方式完成：
SVD直接解決了問題的核心：分解協方差矩陣并用小的正特征值替換負特征值
epsilon。這將使您的矩陣回到 PD 域，從而引入最小的誤差。如果 SVD 的計算成本太高，另一種數值方法就是添加
epsilon * np.eye(D)到協方差矩陣中。epsilon這具有添加到每個特征值的效果。同樣，如果epsilon足夠小，就不會引入太大的錯誤。
如果您采用最后一種方法，您將需要告訴gaussian_kde修改其協方差矩陣。我發現這是一種相對干凈的方法：只需將此類添加到代碼庫中并替換為gaussian_kde（GaussianKde在我這邊測試過，似乎工作正常）。
? ? class GaussianKde(gaussian_kde):
? ? ? ? """
? ? ? ? Drop-in replacement for gaussian_kde that adds the class attribute EPSILON
? ? ? ? to the covmat eigenvalues, to prevent exceptions due to numerical error.
? ? ? ? """
? ??
? ? ? ? EPSILON = 1e-10? # adjust this at will
? ??
? ? ? ? def _compute_covariance(self):
? ? ? ? ? ? """Computes the covariance matrix for each Gaussian kernel using
? ? ? ? ? ? covariance_factor().
? ? ? ? ? ? """
? ? ? ? ? ? self.factor = self.covariance_factor()
? ? ? ? ? ? # Cache covariance and inverse covariance of the data
? ? ? ? ? ? if not hasattr(self, '_data_inv_cov'):
? ? ? ? ? ? ? ? self._data_covariance = np.atleast_2d(np.cov(self.dataset, rowvar=1,
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?bias=False,
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?aweights=self.weights))
? ? ? ? ? ? ? ? # we're going the easy way here
? ? ? ? ? ? ? ? self._data_covariance += self.EPSILON * np.eye(
? ? ? ? ? ? ? ? ? ? len(self._data_covariance))
? ? ? ? ? ? ? ? self._data_inv_cov = np.linalg.inv(self._data_covariance)
? ??
? ? ? ? ? ? self.covariance = self._data_covariance * self.factor**2
? ? ? ? ? ? self.inv_cov = self._data_inv_cov / self.factor**2
? ? ? ? ? ? L = np.linalg.cholesky(self.covariance * 2 * np.pi)
? ? ? ? ? ? self._norm_factor = 2*np.log(np.diag(L)).sum()? # needed for scipy 1.5.2
? ? ? ? ? ? self.log_det = 2*np.log(np.diag(L)).sum()? # changed var name on 1.6.2
解釋
如果您的錯誤類似，但不完全一樣，或者有人感到好奇，這是我遵循的過程，希望它有所幫助。
異常堆棧指定錯誤是在 Cholesky 分解過程中觸發的。就我而言，這是方法內的這一行
_compute_covariance。事實上，在異常發生后，通過檢查
covariance和屬性的特征值表明具有接近于零的負特征值，因此其逆特征值很大。由于 Cholesky 期望 PD 矩陣，這會引發錯誤。inv_covnp.eighcovariance此時，我們可以嚴重懷疑微小的負特征值是一個數值誤差，因為協方差矩陣是 PSD。事實上，當協方差矩陣最初是根據已傳遞給構造函數的數據集計算時，就會出現錯誤源。就我而言，協方差矩陣產生至少 2 個幾乎相同的列，這導致由于數值誤差而產生剩余負特征值。
添加回答
舉報
0/150
提交
取消