建设校园门户网站方案,seo教程排名第一,如何建立一个论坛网站,什么 门户网站Anomaly detection
异常检测是机器学习中比较常见的应用#xff0c;它主要用于非监督学习问题#xff0c;从某些角度看#xff0c; 它又类似于一些监督学习问题。 什么是异常检测#xff1f;来看几个例子#xff1a; 例1. 假设是飞机引擎制造商#xff0c; 要对引擎进行…Anomaly detection
异常检测是机器学习中比较常见的应用它主要用于非监督学习问题从某些角度看 它又类似于一些监督学习问题。 什么是异常检测来看几个例子 例1. 假设是飞机引擎制造商 要对引擎进行质量控制测试现在测量了飞机引擎的一些特征变量比如引擎运转时产生的热量,引擎的震动然后得到了两个特征的无标签数据集 现在有了新的特征变量xtest ,我们希望这个新飞机引擎是否有某种异常。就是异常检测。
在这次练习中实现一个异常检测算法来检测服务器计算机中的异常行为。这些特性测量每个服务器的响应速度 (mb/s)和响应延迟(ms)。 收集到了m307的样本无标签。相信其中绝大多数样本是正常的但还是有一小部分的样本是异常的。
mat loadmat(ex8data1.mat)
print(mat.keys())
X mat[X]
Xval, yval mat[Xval], mat[yval]
X.shape, Xval.shape, yval.shapedef gaussian(X, mu, sigma2):m, n X.shapeif np.ndim(sigma2) 1:sigma2 np.diag(sigma2)norm 1. / (np.power((2 * np.pi), n / 2) * np.sqrt(np.linalg.det(sigma2)))exp np.zeros((m, 1))for row in range(m):xrow X[row]exp[row] np.exp(-0.5 * ((xrow - mu).T).dot(np.linalg.inv(sigma2)).dot(xrow - mu))return norm * expGaussian distribution
高斯分布也称为正态分布变量x 符合高斯分布 xN(µ,σ2)x~N(µ, σ^{2})x N(µ,σ2) 则其概率密度函数为: 利用已有的数据来预测总体中的µ和σ2的计算方法如下:利用已有的数据来预测总体中的µ和σ^{2}的计算方法如下:利用已有的数据来预测总体中的µ和σ2的计算方法如下: 一旦我们获得了平均值和方差的估计值给定新的一个训练实例根据模型计算 p(x) 当p(x) ε, 为异常。
Estimating parameters for a Gaussian
def get_gaussian_params(X, useMultivariate):mu X.mean(axis0)if useMultivariate:sigma2 ((X - mu).T (X - mu)) / len(X)else:sigma2 X.var(axis0, ddof0) # 样本方差return mu, sigma2def plot_contours(mu, sigma2):画出高斯概率分布的图在三维中是一个上凸的曲面。投影到平面上则是一圈圈的等高线。delta .3 # 注意delta不能太小否则会生成太多的数据导致矩阵相乘会出现内存错误。x np.arange(0, 30, delta)y np.arange(0, 30, delta)# 这部分要转化为X形式的坐标矩阵也就是一列是横坐标一列是纵坐标# 然后才能传入gaussian中求解得到每个点的概率值xx, yy np.meshgrid(x, y)points np.c_[xx.ravel(), yy.ravel()] # 按列合并一列横坐标一列纵坐标z gaussian(points, mu, sigma2)z z.reshape(xx.shape) # 这步骤不能忘cont_levels [10 ** h for h in range(-20, 0, 3)]plt.contour(xx, yy, z, cont_levels) # 这个levels是作业里面给的参考,或者通过求解的概率推出来。plt.title(Gaussian Contours, fontsize16)First contours without using multivariate gaussian
plot_data()
useMV False
plot_contours(*get_gaussian_params(X, useMV))# Then contours with multivariate gaussian:
plot_data()
useMV True
# *表示解元组
plot_contours(*get_gaussian_params(X, useMV))Selecting the threshold, ε
确定哪些例子是异常的一种方法是通过一组交叉验证集选择一个好的阈值 ε 。
在这部分的练习中,您将实现一个算法使用交叉验证集的F1分数来选择合理的阈值 ε
tp means true positives是异常值并且我们的模型预测成异常值了即真的异常值。 fp means false positives是正常值但模型把它预测成异常值即假的异常值。 fn means false negatives是异常值但是模型把它预测成正常值即假的正常值。
precision 表示你预测为positive的样本中有多少是真的positive的样本。 recall 表示实际有多少positive的样本而你成功预测出多少positive的样本。
def select_threshold(yval, pval):def computeF1(yval, pval):m len(yval)tp float(len([i for i in range(m) if pval[i] and yval[i]]))fp float(len([i for i in range(m) if pval[i] and not yval[i]]))fn float(len([i for i in range(m) if not pval[i] and yval[i]]))prec tp / (tp fp) if (tp fp) else 0rec tp / (tp fn) if (tp fn) else 0F1 2 * prec * rec / (prec rec) if (prec rec) else 0return F1epsilons np.linspace(min(pval), max(pval), 1000)bestF1, bestEpsilon 0, 0for e in epsilons:pval_ pval ethisF1 computeF1(yval, pval_)if thisF1 bestF1:bestF1 thisF1bestEpsilon ereturn bestF1, bestEpsilonmu, sigma2 get_gaussian_params(X, useMultivariateFalse)
pval gaussian(Xval, mu, sigma2)
bestF1, bestEpsilon select_threshold(yval, pval)
# (0.8750000000000001, 8.999852631901397e-05)y gaussian(X, mu, sigma2) # X的概率
xx np.array([X[i] for i in range(len(y)) if y[i] bestEpsilon])
xx # 离群点plot_data()
plot_contours(mu, sigma2)
plt.scatter(xx[:, 0], xx[:, 1], s80, facecolorsnone, edgecolorsr)# High dimensional dataset
mat loadmat(ex8data2.mat)
X2 mat[X]
Xval2, yval2 mat[Xval], mat[yval]
X2.shape # (1000, 11)mu, sigma2 get_gaussian_params(X2, useMultivariateFalse)
ypred gaussian(X2, mu, sigma2)yval2pred gaussian(Xval2, mu, sigma2)# You should see a value epsilon of about 1.38e-18, and 117 anomalies found.
bestF1, bestEps select_threshold(yval2, yval2pred)
anoms [X2[i] for i in range(X2.shape[0]) if ypred[i] bestEps]
bestEps, len(anoms)Recommender Systems
Movie ratings dataset
mat loadmat(ex8_movies.mat)
print(mat.keys())
Y, R mat[Y], mat[R]
nm, nu Y.shape # Y中0代表用户没有评分
nf 100
Y.shape, R.shape
# ((1682, 943), (1682, 943))Y[0].sum() / R[0].sum() # 分子代表第一个电影的总分数分母代表这部电影有多少评分数据# Visualize the ratings matrix
fig plt.figure(figsize(8,8*(1682./943.)))
plt.imshow(Y, cmaprainbow)
plt.colorbar()
plt.ylabel(Movies (%d) % nm, fontsize20)mat loadmat(ex8_movieParams.mat)
X mat[X]
Theta mat[Theta]
nu int(mat[num_users])
nm int(mat[num_movies])
nf int(mat[num_features])
# For now, reduce the data set size so that this runs faster
nu 4
nm 5
nf 3
X X[:nm, :nf]
Theta Theta[:nu, :nf]
Y Y[:nm, :nu]
R R[:nm, :nu]Collaborative filtering learning algorithm
# 展开参数
def serialize(X, Theta):return np.r_[X.flatten(),Theta.flatten()]# 提取参数
def deserialize(seq, nm, nu, nf):return seq[:nm*nf].reshape(nm, nf), seq[nm*nf:].reshape(nu, nf)Collaborative filtering cost function
def co_fi_cost_func(params, Y, R, nm, nu, nf, l0):params : 拉成一维之后的参数向量(X, Theta)Y : 评分矩阵 (nm, nu)R 0-1矩阵表示用户对某一电影有无评分nu : 用户数量nm : 电影数量nf : 自定义的特征的维度l : lambda for regularizationX, Theta deserialize(params, nm, nu, nf)# (XTheta)*R含义如下 因为XTheta是我们用自定义参数算的评分但是有些电影本来是没有人# 评分的存储在R中0-1表示。将这两个相乘得到的值就是我们要的已经被评分过的电影的预测分数。error 0.5 * np.square((X Theta.T - Y) * R).sum()reg1 .5 * l * np.square(Theta).sum()reg2 .5 * l * np.square(X).sum()return error reg1 reg2co_fi_cost_func(serialize(X,Theta),Y,R,nm,nu,nf),co_fi_cost_func(serialize(X,Theta),Y,R,nm,nu,nf,1.5)
Collaborative filtering gradient
def check_gradient(params, Y, myR, nm, nu, nf, l0.):Lets check my gradient computation real quickprint(Numerical Gradient \t cofiGrad \t\t Difference)# 分析出来的梯度grad co_fi_gradient(params, Y, myR, nm, nu, nf, l)# 用 微小的e 来计算数值梯度。e 0.0001nparams len(params)e_vec np.zeros(nparams)# Choose 10 random elements of param vector and compute the numerical gradient# 每次只能改变e_vec中的一个值并在计算完数值梯度后要还原。for i in range(10):idx np.random.randint(0, nparams)e_vec[idx] eloss1 co_fi_cost_func(params - e_vec, Y, myR, nm, nu, nf, l)loss2 co_fi_cost_func(params e_vec, Y, myR, nm, nu, nf, l)numgrad (loss2 - loss1) / (2 * e)e_vec[idx] 0diff np.linalg.norm(numgrad - grad[idx]) / np.linalg.norm(numgrad grad[idx])print(%0.15f \t %0.15f \t %0.15f % (numgrad, grad[idx], diff))print(Checking gradient with lambda 0...)
check_gradient(serialize(X, Theta), Y, R, nm, nu, nf)
print(\nChecking gradient with lambda 1.5...)
check_gradient(serialize(X, Theta), Y, R, nm, nu, nf, l1.5)Learning movie recommendations
movies [] # 包含所有电影的列表
with open(data/movie_ids.txt, r, encodingutf-8) as f:for line in f:
# movies.append( .join(line.strip().split( )[1:]))movies.append( .join(line.strip().split( )[1:]))my_ratings np.zeros((1682, 1))my_ratings[0] 4
my_ratings[97] 2
my_ratings[6] 3
my_ratings[11] 5
my_ratings[53] 4
my_ratings[63] 5
my_ratings[65] 3
my_ratings[68] 5
my_ratings[182] 4
my_ratings[225] 5
my_ratings[354] 5for i in range(len(my_ratings)):if my_ratings[i] 0:print(my_ratings[i], movies[i])mat loadmat(data/ex8_movies.mat)
Y, R mat[Y], mat[R]
Y.shape, R.shapeY np.c_[Y, my_ratings] # (1682, 944)
R np.c_[R, my_ratings!0] # (1682, 944)
nm, nu Y.shape
nf 10 # 我们使用10维的特征向量def normalizeRatings(Y, R):The mean is only counting movies that were ratedYmean (Y.sum(axis1) / R.sum(axis1)).reshape(-1,1)
# Ynorm (Y - Ymean)*R # 这里也要注意不要归一化未评分的数据Ynorm (Y - Ymean)*R # 这里也要注意不要归一化未评分的数据return Ynorm, YmeanYnorm, Ymean normalizeRatings(Y, R)
Ynorm.shape, Ymean.shape
# ((1682, 944), (1682, 1))X np.random.random((nm, nf))
Theta np.random.random((nu, nf))
params serialize(X, Theta)
l 10Recommendations
import scipy.optimize as opt
res opt.minimize(funco_fi_cost_func(),x0params,args(Y, R, nm, nu, nf, l),methodTNC,jacco_fi_gradient(),options{maxiter: 100})ret res.xfit_X, fit_Theta deserialize(ret, nm, nu, nf)# Recommendations
# 所有用户的剧场分数矩阵
pred_mat fit_X fit_Theta.T# 最后一个用户的预测分数 也就是我们刚才添加的用户
pred pred_mat[:,-1] Ymean.flatten()pred_sorted_idx np.argsort(pred)[::-1] # 排序并翻转使之从大到小排列print(Top recommendations for you:)
for i in range(10):print(Predicting rating %0.1f for movie %s. \% (pred[pred_sorted_idx[i]], movies[pred_sorted_idx[i]]))print(\nOriginal ratings provided:)
for i in range(len(my_ratings)):if my_ratings[i] 0:print(Rated %d for movie %s. % (my_ratings[i], movies[i]))
