我有这个有异常值的数据。我怎样才能找到马氏距离并用它来消除异常值。
问问题
11541 次
3 回答
44
首先让我提出一些一般准则：
- 实际上，如果你有很多特征和更少的样本，马氏算法往往会给出误导性的结果（你可以自己尝试一下），所以你拥有的特征越多，你应该提供的样本越多。
- 协方差矩阵必须是对称的和正定的才能使算法有效，因此您应该在继续之前进行检查。
正如已经提到的，欧几里得度量未能找到正确的距离，因为它试图获得普通的直线距离。因此，如果我们有变量的多维空间，两个点可能看起来与Mean的距离相同，但其中一个点距离数据云很远（即它是一个异常值）。
解决方案是马氏距离，它通过采用变量的特征向量而不是原始轴来实现类似于特征缩放的效果。
它应用以下公式：
在哪里：
x是观察到的距离；m是观测值的平均值；S是协方差矩阵。
复习：
协方差表示两个变量之间关系的方向（即正、负或零），因此它显示了一个变量如何与其他变量的变化相关的强度。
执行
考虑这个6x3数据集，其中每一行代表一个样本，每一列代表给定样本的一个特征：
首先，我们需要为每个样本的特征创建一个协方差矩阵，这就是为什么我们在numpy.covrowvar函数中将参数设置为，所以现在每列代表一个变量：False
covariance_matrix = np.cov(data, rowvar=False)
# data here looks similar to the above table
# in the picture
接下来，我们找到协方差矩阵的逆：
inv_covariance_matrix = np.linalg.inv(covariance_matrix)
但在继续之前，如上所述，我们应该检查矩阵及其逆矩阵是否是对称的和正定的。我们使用这个Cholesky 分解算法，幸运的是，它已经在numpy.linalg.cholesky中实现了：
def is_pos_def(A):
if np.allclose(A, A.T):
try:
np.linalg.cholesky(A)
return True
except np.linalg.LinAlgError:
return False
else:
return False
然后，我们找到m每个特征上的变量的平均值（我应该说维度）并将它们保存在一个数组中，如下所示：
vars_mean = []
for i in range(data.shape[0]):
vars_mean.append(list(data.mean(axis=0)))
# axis=0 means each column in the 2D array
请注意，我重复每一行只是为了利用矩阵减法，如下所示。
接下来，我们找到x - m（即微分），但由于我们已经有了矢量化vars_mean，我们需要做的就是：
diff = data - vars_mean
# here we subtract the mean of feature
# from each feature of each example
最后，应用这样的公式：
md = []
for i in range(len(diff)):
md.append(np.sqrt(diff[i].dot(inv_covariance_matrix).dot(diff[i])))
请注意以下事项：
- 协方差矩阵逆的维数为：
number_of_features x number_of_features - 矩阵的维度
diff类似于原始数据矩阵：number_of_examples x number_of_features - 因此，每个
diff[i]（即行）是1 x number_of_features。 - 所以根据矩阵乘法规则，结果矩阵
diff[i].dot(inv_covariance_matrix)将是1 x number_of_features; 当我们再次乘以diff[i];numpy自动将后者视为列矩阵（即number_of_features x 1）；所以最终的结果会变成一个单一的值（即不需要转置）。
为了检测异常值，我们应该指定一个threshold; 但是由于 Mahalanobis 距离的平方遵循卡方分布，其自由度 = 数据集中的特征数，所以我们可以选择 0.1 的阈值，然后我们可以使用chi2.cdf方法Scipy，如下所示：
1 - chi2.cdf(square_of_mahalanobis_distances, degree_of_freedom)
因此，任何具有 (1 - 卡方 CDF) 小于或等于阈值的点都可以归类为异常值。
放在一起
import numpy as np
def create_data(examples=50, features=5, upper_bound=10, outliers_fraction=0.1, extreme=False):
'''
This method for testing (i.e. to generate a 2D array of data)
'''
data = []
magnitude = 4 if extreme else 3
for i in range(examples):
if (examples - i) <= round((float(examples) * outliers_fraction)):
data.append(np.random.poisson(upper_bound ** magnitude, features).tolist())
else:
data.append(np.random.poisson(upper_bound, features).tolist())
return np.array(data)
def MahalanobisDist(data, verbose=False):
covariance_matrix = np.cov(data, rowvar=False)
if is_pos_def(covariance_matrix):
inv_covariance_matrix = np.linalg.inv(covariance_matrix)
if is_pos_def(inv_covariance_matrix):
vars_mean = []
for i in range(data.shape[0]):
vars_mean.append(list(data.mean(axis=0)))
diff = data - vars_mean
md = []
for i in range(len(diff)):
md.append(np.sqrt(diff[i].dot(inv_covariance_matrix).dot(diff[i])))
if verbose:
print("Covariance Matrix:\n {}\n".format(covariance_matrix))
print("Inverse of Covariance Matrix:\n {}\n".format(inv_covariance_matrix))
print("Variables Mean Vector:\n {}\n".format(vars_mean))
print("Variables - Variables Mean Vector:\n {}\n".format(diff))
print("Mahalanobis Distance:\n {}\n".format(md))
return md
else:
print("Error: Inverse of Covariance Matrix is not positive definite!")
else:
print("Error: Covariance Matrix is not positive definite!")
def is_pos_def(A):
if np.allclose(A, A.T):
try:
np.linalg.cholesky(A)
return True
except np.linalg.LinAlgError:
return False
else:
return False
data = create_data(15, 3, 10, 0.1)
print("data:\n {}\n".format(data))
MahalanobisDist(data, verbose=True)
结果
data:
[[ 12 7 9]
[ 9 16 7]
[ 14 11 10]
[ 14 5 5]
[ 12 8 7]
[ 8 8 10]
[ 9 14 8]
[ 12 12 10]
[ 18 10 6]
[ 6 12 11]
[ 4 12 15]
[ 5 13 10]
[ 8 9 8]
[106 116 97]
[ 90 116 114]]
Covariance Matrix:
[[ 980.17142857 1143.62857143 1035.6 ]
[1143.62857143 1385.11428571 1263.12857143]
[1035.6 1263.12857143 1170.74285714]]
Inverse of Covariance Matrix:
[[ 0.03021777 -0.03563241 0.0117146 ]
[-0.03563241 0.08684092 -0.06217448]
[ 0.0117146 -0.06217448 0.05757261]]
Variables Mean Vector:
[[21.8, 24.6, 21.8], [21.8, 24.6, 21.8], [21.8, 24.6, 21.8], [21.8, 24.6, 21.8], [21.8, 24.6, 21.8], [21.8, 24.6, 21.8], [21.8, 24.6, 21.8], [21.8, 24.6, 21.8], [21.8, 24.6, 21.8], [21.8, 24.6, 21.8], [21.8, 24.6, 21.8], [21.8, 24.6, 21.8], [21.8, 24.6, 21.8], [21.8, 24.6, 21.8], [21.8, 24.6, 21.8]]
Variables - Variables Mean Vector:
[[ -9.8 -17.6 -12.8]
[-12.8 -8.6 -14.8]
[ -7.8 -13.6 -11.8]
[ -7.8 -19.6 -16.8]
[ -9.8 -16.6 -14.8]
[-13.8 -16.6 -11.8]
[-12.8 -10.6 -13.8]
[ -9.8 -12.6 -11.8]
[ -3.8 -14.6 -15.8]
[-15.8 -12.6 -10.8]
[-17.8 -12.6 -6.8]
[-16.8 -11.6 -11.8]
[-13.8 -15.6 -13.8]
[ 84.2 91.4 75.2]
[ 68.2 91.4 92.2]]
Mahalanobis Distance:
[1.3669401667524865, 2.1796331318432967, 0.7470525416547134, 1.6364973119931507, 0.8351423113609481, 0.9128858131134882, 1.397144258271586, 0.35603382066414996, 1.4449501739129382, 0.9668775289588046, 1.490503433100514, 1.4021488309805878, 0.4500345257064412, 3.239353067840299, 3.260149280200771]
于 2018-08-17T11:23:46.543 回答
5
在多元数据中，如果变量之间存在协方差（即在您的情况下为 X、Y、Z），欧几里得距离将失败。
因此，马氏距离所做的是，
它将变量转换为不相关的空间。
使每个变量的方差等于 1。
然后计算简单的欧几里得距离。
我们可以如下计算每个数据样本的马氏距离，
在这里，我提供了python代码并添加了注释，以便您理解代码。
import numpy as np
data= np.matrix([[1, 2, 3, 4, 5, 6, 7, 8],[1, 4, 9, 16, 25, 36, 49, 64],[1, 4, 9, 16, 25, 16, 49, 64]])
def MahalanobisDist(data):
covariance_xyz = np.cov(data) # calculate the covarince matrix
inv_covariance_xyz = np.linalg.inv(covariance_xyz) #take the inverse of the covarince matrix
xyz_mean = np.mean(data[0]),np.mean(data[1]),np.mean(data[2])
x_diff = np.array([x_i - xyz_mean[0] for x_i in x]) # take the diffrence between the mean of X variable the sample
y_diff = np.array([y_i - xyz_mean[1] for y_i in y]) # take the diffrence between the mean of Y variable the sample
z_diff = np.array([z_i - xyz_mean[2] for z_i in z]) # take the diffrence between the mean of Z variable the sample
diff_xyz = np.transpose([x_diff, y_diff, z_diff])
md = []
for i in range(len(diff_xyz)):
md.append(np.sqrt(np.dot(np.dot(np.transpose(diff_xyz[i]),inv_covariance_xyz),diff_xyz[i]))) #calculate the Mahalanobis Distance for each data sample
return md
def MD_removeOutliers(data):
MD = MahalanobisDist(data)
threshold = np.mean(MD) * 1.5 # adjust 1.5 accordingly
outliers = []
for i in range(len(MD)):
if MD[i] > threshold:
outliers.append(i) # index of the outlier
return np.array(outliers)
print(MD_removeOutliers(data))
希望这可以帮助。
参考，
于 2017-10-19T15:34:11.760 回答
0
正如@Yahya 之前指出的那样，为了计算 MD，cov .matrix 应该是正的和半定的。这些是取 cov 的倒数的必要条件。矩阵 。矩阵的逆使用矩阵的倒数determinant。您可能无法计算矩阵的行列式这一事实可能表明数据集本身存在更深层次的问题。这可能是因为数据集中的两列或多列是相关的。您最好在尝试计算 MD 之前删除这些对。使用Pseudo Inverse可以是另一种选择。
于 2022-01-11T03:20:07.250 回答