机器学习实践七----异常检测和推荐系统

最新推荐文章于 2024-04-01 00:44:51 发布

破岩

最新推荐文章于 2024-04-01 00:44:51 发布

阅读量492

点赞数

分类专栏：机器学习文章标签：人工智能机器学习

本文链接：https://blog.csdn.net/stevenhenry/article/details/104622379

版权

机器学习专栏收录该内容

8 篇文章 0 订阅

订阅专栏

Anomaly detection

异常检测是机器学习中比较常见的应用，它主要用于非监督学习问题，从某些角度看，它又类似于一些监督学习问题。
什么是异常检测？来看几个例子：
例1. 假设是飞机引擎制造商，要对引擎进行质量控制测试，现在测量了飞机引擎的一些特征变量，比如引擎运转时产生的热量,引擎的震动，然后得到了两个特征的无标签数据集，现在有了新的特征变量xtest ,我们希望这个新飞机引擎是否有某种异常。就是异常检测。

在这次练习中，实现一个异常检测算法来检测服务器计算机中的异常行为。这些特性测量每个服务器的响应速度 (mb/s)和响应延迟(ms)。
收集到了m=307的样本，无标签。相信其中绝大多数样本是正常的，但还是有一小部分的样本是异常的。

mat = loadmat('ex8data1.mat')
print(mat.keys())
X = mat['X']
Xval, yval = mat['Xval'], mat['yval']
X.shape, Xval.shape, yval.shape

def gaussian(X, mu, sigma2):
    m, n = X.shape
    if 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 * exp

Gaussian distribution

高斯分布，也称为正态分布，变量x 符合高斯分布 $x~N(µ, σ^{2})$ 则其概率密度函数为:
在这里插入图片描述
$利用已有的数据来预测总体中的µ和σ^{2}的计算方法如下:$

一旦我们获得了平均值和方差的估计值，给定新的一个训练实例，根据模型计算 p(x)：

当p(x) < ε, 为异常。

Estimating parameters for a Gaussian

def get_gaussian_params(X, useMultivariate):
    mu = X.mean(axis=0)
    if useMultivariate:
        sigma2 = ((X - mu).T @ (X - mu)) / len(X)
    else:
        sigma2 = X.var(axis=0, ddof=0)  # 样本方差

    return mu, sigma2

def 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', fontsize=16)

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 0
        rec = tp / (tp + fn) if (tp + fn) else 0
        F1 = 2 * prec * rec / (prec + rec) if (prec + rec) else 0
        return F1

    epsilons = np.linspace(min(pval), max(pval), 1000)
    bestF1, bestEpsilon = 0, 0
    for e in epsilons:
        pval_ = pval < e
        thisF1 = computeF1(yval, pval_)
        if thisF1 > bestF1:
            bestF1 = thisF1
            bestEpsilon = e

    return bestF1, bestEpsilon

mu, sigma2 = get_gaussian_params(X, useMultivariate=False)
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], s=80, facecolors='none', edgecolors='r')

# 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, useMultivariate=False)
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, cmap='rainbow')
plt.colorbar()
plt.ylabel('Movies (%d) ' % nm, fontsize=20)

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, l=0):
    """
    params : 拉成一维之后的参数向量(X, Theta)
    Y : 评分矩阵 (nm, nu)
    R ：0-1矩阵，表示用户对某一电影有无评分
    nu : 用户数量
    nm : 电影数量
    nf : 自定义的特征的维度
    l : lambda for regularization
    """
    X, Theta = deserialize(params, nm, nu, nf)

    # (X@Theta)*R含义如下： 因为X@Theta是我们用自定义参数算的评分，但是有些电影本来是没有人
    # 评分的，存储在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 + reg2


co_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, l=0.):
    """
    Let's check my gradient computation real quick
    """
    print('Numerical Gradient \t cofiGrad \t\t Difference')

    # 分析出来的梯度
    grad = co_fi_gradient(params, Y, myR, nm, nu, nf, l)

    # 用 微小的e 来计算数值梯度。
    e = 0.0001
    nparams = 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] = e
        loss1 = 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] = 0
        diff = 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, l=1.5)

Learning movie recommendations

movies = []  # 包含所有电影的列表
with open('data/movie_ids.txt', 'r', encoding='utf-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] = 5

for 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.shape

Y = 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 rated
    """
    Ymean = (Y.sum(axis=1) / R.sum(axis=1)).reshape(-1,1)
#     Ynorm = (Y - Ymean)*R  # 这里也要注意不要归一化未评分的数据
    Ynorm = (Y - Ymean)*R  # 这里也要注意不要归一化未评分的数据
    return Ynorm, Ymean

Ynorm, 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 = 10

Recommendations

import scipy.optimize as opt
res = opt.minimize(fun=co_fi_cost_func(),
                   x0=params,
                   args=(Y, R, nm, nu, nf, l),
                   method='TNC',
                   jac=co_fi_gradient(),
                   options={'maxiter': 100})

ret = res.x

fit_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]))

破岩

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
机器学习实践七----异常检测和推荐系统

Anomaly detectionGaussian distributionEstimating parameters for a GaussianSelecting the threshold, εHigh dimensional datasetRecommender SystemsMovie ratings datasetCollaborative filtering learn...
复制链接

扫一扫

专栏目录