书接上文,继续讨论基于多元正态分布的异常检测算法。
现在有一个包含了m个数据的训练集,其中的每个样本都是一个n维数据:
可以通过下面的函数判断一个样本是否是异常的:
我们的目的是设法根据训练集求得μ和σ,以得到一个确定的多元分正态布模型。具体来说,通过最大似然估计量可以得出下面的结论:
其中Σ是协方差对角矩阵,最终求得的多元正态分布模型可以写成:
关于最大似然估计量、协方差矩阵和多元正态分布最大似然估计的具体推导过程,可参考:
代码:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
def create_data(model='train', count=200):
'''
构造2维训练集
:param model: train:训练集, test:测试集
:param count: 样本数量
:return: X1:第1纬度的列, X2:第2维度的列表
'''
np.random.seed(21) # 设置seed使每次生成的随机数都相等
X1 = np.random.normal(1.7, 0.036, count) # 生成200个符合正态分布的身高数据
low, high = -0.01, 0.01
if model == 'test':
low, high = -0.05, 0.05
# 设置身高对应的鞋码,正常鞋码=身高/7±0.01
X2 = X1 / 7 + np.random.uniform(low=low, high=high, size=len(X1))
return X1, X2
def plot_train(X1, X2):
'''
可视化训练集
:param X1: 训练集的第1维度
:param X2: 训练集的第2维度
'''
fig = plt.figure(figsize=(10, 4))
plt.subplots_adjust(hspace=0.5) # 调整子图之间的上下边距
# 数据的散点图
fig.add_subplot(2, 1, 1)
plt.scatter(X1, X2, color='k', s=3., label='训练数据')
plt.legend(loc='upper left')
plt.xlabel('身高')
plt.ylabel('脚长')
plt.title('数据分布')
# 身高维度的直方图
fig.add_subplot(2, 2, 3)
plt.hist(X1, bins=40)
plt.xlabel('身高')
plt.ylabel('频度')
plt.title('身高直方图')
# 脚长维度的直方图
fig.add_subplot(2, 2, 4)
plt.hist(X2, bins=40)
plt.xlabel('脚长')
plt.ylabel('频度')
plt.title('脚长直方图')
plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False # 解决中文下的坐标轴负号显示问题
plt.show()
def fit(X_train):
'''
拟合数据,训练模型
:param X_train: 训练集
:return: mu:均值, sigma:方差
'''
global mu, sigma
X = np.mat(X_train.T)
m, n = X.shape
mu = np.mean(X, axis=1) # 计算均值μ,axis=1表示对每一个子数组计算均值
sigma = np.mat(np.cov(X)) # 计算Σ,等同于(X - mu) * (X - mu).T / len(X_train)
def gaussian(X_test):
'''
计算正态分布
:param X_test: 测试集
:return: 数据集的密度值
'''
global mu, sigma
m, n = np.shape(X_test)
sig_det = np.linalg.det(sigma) # 计算det(Σ)
sig_inv = np.linalg.inv(sigma) # Σ的逆矩阵
r = []
for x in X_test:
x = np.mat(x).T - mu
g = np.exp(-x.T * sig_inv * x / 2) * ((2 * np.pi) ** (-n / 2) * (sig_det ** (-0.5)))
r.append(g[0, 0])
return r
def sel_epsilon(X_train):
'''
选择合适的ε
:param X_train:
:return: ε
'''
g_val = gaussian(X_train)
return np.min(g_val) + 0.0001
def predict(X):
'''
检测训练集中的数据是否是正常数据
:param X: 待预测的数据
:return: P1:数据的密度值, P2:数据的异常检测结果,True正常,False异常
'''
P1 = gaussian(X) # 数据的密度值
P2 = [p > epsilon for p in P1] # 数据的异常检测结果,True正常,False异常
return P1, P2
def plot_predict(X):
'''可视化异常检测结果 '''
P1, P2 = predict(X)
normals_idx = [i for i, t in enumerate(P2) if t == True] # 正常数据的索引
outliers_idx = [i for i, t in enumerate(P2) if t == False] # 异常数据的索引
normals_x = np.array([X[i] for i in normals_idx]) # 正常数据
outliers_x = np.array([X[i] for i in outliers_idx]) # 异常数据
fig1 = plt.figure(num='fig1') # 散点图
ax = Axes3D(fig1)
ax.scatter(normals_x[:,0], normals_x[:,1],
[P1[i] for i in normals_idx], label='正常数据')
ax.scatter(outliers_x[:,0], outliers_x[:,1],
[P1[i] for i in outliers_idx], c='r', marker='^', label='异常数据')
ax.set_title('共有{}个异常数据'.format(len(outliers_idx)))
ax.axis('tight') # 让坐标轴的比例尺适应数据量
ax.set_xlabel('身高')
ax.set_ylabel('脚长')
ax.set_zlabel('p(x)')
ax.legend(loc='upper left')
n = 100
xs, ys = np.meshgrid(np.linspace(min(X1_train), max(X1_train), n),
np.linspace(min(X2_train), max(X2_train), n))
zs = [gaussian(np.c_[xs[i], ys[i]]) for i in range(n)]
fig2 = plt.figure(num='fig2')
ax = Axes3D(fig2)
# 3维空间的拟合曲面
ax.plot_surface(xs, ys, zs, alpha=0.5, cmap=plt.get_cmap('rainbow'))
ax.scatter(normals_x[:, 0], normals_x[:, 1],
[P1[i] for i in normals_idx], label='正常数据')
ax.scatter(outliers_x[:, 0], outliers_x[:, 1],
[P1[i] for i in outliers_idx], c='r', marker='^', label='异常数据')
ax.axis('tight') # 让坐标轴的比例尺适应数据量
ax.set_xlabel('身高')
ax.set_ylabel('脚长')
ax.set_zlabel('p(x)')
ax.legend(loc='upper left')
fig3 = plt.figure(num='fig3')
plt.scatter(normals_x[:, 0], normals_x[:, 1], s=30., c='k', label='正常数据')
plt.scatter(outliers_x[:, 0], outliers_x[:, 1], c='r', marker='^', s=30., label='异常数据')
plt.contour(xs, ys, zs, 80, alpha=0.5) # 等高线
plt.axis('tight') # 让坐标轴的比例尺适应数据量
plt.xlabel('身高')
plt.ylabel('脚长')
plt.legend(loc='upper left')
plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False # 解决中文下的坐标轴负号显示问题
plt.show()
if __name__ == '__main__':
X1_train, X2_train = create_data()
plot_train(X1_train, X2_train)
X_train = np.c_[X1_train, X2_train]
mu, sigma = np.mat([]), np.mat([])
fit(X_train)
epsilon = sel_epsilon(X_train)
X_test = np.c_[create_data(model='test', count=20)]
X = np.r_[X_train, X_test]
plot_predict(X)
为了简单起见,我们认为X_train 中的数据全部是正常数据。fit(X_train)计算多元正态分布的模型参数,gaussian(X_test)根据目标函数计算样本的多元正态分布密度值。在知道了算法原理的请看下,fit(X_train)和gaussian(X_test)都毫无神秘性可言。
predict(X)将对X中的所有样本进行检测,并返回X对应的检测结果列表,列表中的元素是一个二元组,第一个元素记录x(i)是否是正常数据,第二个元素记录p(x(i);μ,Σ)。由于已经假设了X_train中全部是正常数据,因此这里选择X_train中最小的密度值作为ε。
X_test中的20个测试数据是可能的异常样本。plot_predict(X)展示了空间样本点、空间拟合曲面和等高线:
作者:我是8位的
出处:http://www.cnblogs.com/bigmonkey
本文以学习、研究和分享为主,如需转载,请联系本人,标明作者和出处,非商业用途!
扫描二维码关注公作者众号“我是8位的”