task_3异常检测——线性相关方法

2、数据可视化

以breast-cancer-unsupervised-ad数据集为例做一些简单的数据可视化。

#coding:utf-8
#导入warnings包，利用过滤器来实现忽略警告语句。?这个不知道是干什么用的
import warnings
warnings.filterwarnings('ignore')
import numpy as np
import pandas as pd 
import matplotlib.pyplot as plt
import seaborn as sns
import missingno as msno

# 载入训练集和测试集
f = open('dataverse_files/breast-cancer-unsupervised-ad.csv')
Train_data = pd.read_csv(f)
Train_data.head()
Train_data.tail()
Train_data.describe()# 均值方差数量级别都相差很大
#Train_data.info()#最后一列是标签，其他都是浮点型
numeric_features = ['f' + str(i) for i in range(30)]
numeric_features
# 相关性分析
numeric = Train_data[numeric_features]
numeric#不加标签的数据
correlation  = numeric.corr()
correlation
f,ax = plt.subplots(figsize=(14,14))
sns.heatmap(correlation,square=True)
plt.title('Correlation of Numeric Features with Price',y=1,size=8)
plt.show()

## 3) 每个数字特征的分布可视化
f = pd.melt(Train_data,value_vars=numeric_features)#宽表转换为长表，即行转列
f.head()
#FacetGrid当想要在数据集的子集中分别可视化变量的分布或多个变量之间的关系时，该类非常有用。
#data : DataFrame
#处理后的（“长格式”）dataframe数据，其中每一列都是一个变量（特征），每一行都是一个样本
#row, col, hue : strings 
#定义数据子集的变量，这些变量将在网格的不同方面绘制。请参阅下面*_order参数以控制该变量的级别顺序
#col_wrap : int, optional
#这个意思是图网格列维度限制，比如col_wrap =3，那么在这个画布里最多只能画3列。行不限制，这样就限制了列的个数。
#share{x,y} : bool, ‘col’, or ‘row’ optional
#是否共享x轴或者y轴，就是说如果为真，就共享同一个轴，否则就不共享，默认是都共享，即都为True
#基本工作流程是FacetGrid使用数据集和用于构造网格的变量初始化对象。然后，可以通过调用FacetGrid.map()或将一个或多个绘图函数应用于每个子集 FacetGrid.map_dataframe()。
#最后，可以使用其他方法调整绘图，以执行更改轴标签，使用不同刻度或添加图例等操作
g = sns.FacetGrid(f,col='variable',col_wrap = 6,sharex=False,sharey=False)
g.map(sns.distplot,"value", hist=False, rug=True)#??

#变量两两之间的相关性
sns.set()#在绘制图形之前调用Seaborn的set()函数就可以直接使用其设定好的默认主题进行美化
sns.pairplot(numeric,size=2,kind='scatter',diag_kind='kde')
#plt.savefig('correlation.png')
plt.show()#如果运行则显示一直在忙，不知道是什么情况

#数据降维可视化
from sklearn.manifold import TSNE#t-distributed Stochastic Neighbor Embedding，t分布的随机相邻嵌入。
#具体可参考https://rosefun.blog.csdn.net/article/details/78941241?utm_medium=distribute.pc_relevant.none-task-blog-OPENSEARCH-2.control&depth_1-utm_source=distribute.pc_relevant.none-task-blog-OPENSEARCH-2.control
tsne = TSNE(n_components=2, init='pca', random_state=0)#目的是为了说明，要用pca的方法进行将行降维，
#参数设置具体可参考https://blog.csdn.net/weixin_44530236/article/details/89309046?ops_request_misc=%25257B%252522request%25255Fid%252522%25253A%252522161088776116780274138962%252522%25252C%252522scm%252522%25253A%25252220140713.130102334..%252522%25257D&request_id=161088776116780274138962&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~all~sobaiduend~default-1-89309046.first_rank_v2_pc_rank_v29&utm_term=sklearn.manifold%20import%20TSNE
result = tsne.fit_transform(numeric)#把数据放进去降维
x_min, x_max = np.min(result, 0), np.max(result, 0)
result = (result - x_min) / (x_max - x_min)#归一化处理
label = Train_data['label']
fig = plt.figure(figsize = (7, 7))
#f , ax = plt.subplots()
color = {'o':0, 'n':7}#?这个标记为什么这么设计
for i in range(result.shape[0]):
    plt.text(result[i, 0], result[i, 1], str(label[i]),#plt.text()给图形添加数据标签，通过for循环找到每一个x、y值的相应坐标——result[i,0],result[i,1]，str(label[i]就是要在每个坐标上标注的文字内容
                color=plt.cm.Set1(color[label[i]] / 10.),
                fontdict={'weight': 'bold', 'size': 9})#fontsize则是文字大小
plt.xticks([])#xticks与yticks（设置坐标轴刻度），xlabel与ylabel（设置坐标轴标签）、title（标题）、legend（图例）、xlim与ylim（设置坐标轴数据范围）、
　#　grid（设置网格线）等命令来装饰图形
plt.yticks([])
plt.title('Visualization of data dimension reduction')

3、线性回归

在线性回归中，我们假设不同维度的变量具有一定的相关性，并可以通过一个相关系数矩阵进行衡量。因此对于特定的观测值，可以通过线性方程组来建模。在实际应用中，观测值的数量往往远大于数据的维度，也就是方程个数远远大于未知数的个数，导致线性方程组是一个超定方程，不能直接求解。

• 线性回归通过一系列自变量去预测一个特殊因变量的值，在这种情况下，异常值是根据其他自变量对因变量的影响来定义的，而自变量之间相互关系中的异常则不那么重要。 这里的异常点检测主要用于数据降噪，避免异常点的出现对模型性能的影响，因而这里关注的兴趣点主要是正常值(n)。
• 而我们通常所说的异常检测中并不会对任何变量给与特殊对待，异常值的定义是基于基础数据点的整体分布，这里我们关注的兴趣点主要是异常值(o)。
  广义的回归建模二者都可用

3.1 基于自变量与因变量的线性回归

首先介绍一般的最小二乘方法来求解回归模型

3.1.1 最小二乘法

以一元线性回归为例:
$$Y=\sum _{i=1}^d{a}_i\cdot X_i+a_{d+1}$$
变量Y为因变量，也就是我们要预测的值(一般也就是标签值)；$X_1...X_d$为一系列自变量(样本值)，也就是输入值。系数$a_1...a_{d+1}$为要学习的参数。假设数据共包含$N$个样本，第$j$个样本包含的数据为$x_{j1}...x_{jd}$和$y_j$，带入式(1)如下式所示：

$$y_j=\sum _{i=1}^d{a}_i\cdot x_{ji}+a_{d+1}+{ϵ}_j$$
这里${ϵ}_j$为第$j$个样本的误差。以$Y$ 代表 $N\times 1$ 的因变量矩阵${(y_1...y_N)}^T$，即样本中的真实值；以$U$代表$N\times (d+1)$的自变量矩阵，其中第$j$行为$(x_{j1}...x_{jd},1)$；以$A$ 代表 $(d+1)\times 1$ 的系数矩阵$(a_1...a_{d+1}{)}^T$。则模型可表示为： $$Y=U\cdot A+ϵ$$
定义目标函数为：

$$L(A)=\frac{1}{2}\left||Y-U\cdot A\right|{|}^2$$
目标函数是关于$A$的凸函数，其对$A$求偏导为：

$$\frac{\partial L(A)}{\partial A}=\frac{1}{2}\frac{\partial {\left||Y-U\cdot A\right||}^2}{\partial A}=-U^T(Y-U\cdot A)$$

令$\frac{\partial L(A)}{\partial A}=0$，得到最优参数为：

$$A={\left(U^T\cdot U\right)}^{-1}\cdot \left(U^T\cdot Y\right)$$
这种求解线性回归参数的方法也叫最小二乘法。
最小二乘法要求矩阵 $U^T\cdot U$ 可逆，即$U^T\cdot U$是满秩的。当$U^T\cdot U$不可逆时可以通过两种方法进行参数估计，一种先使用主成分分析等方法来预处理数据，消除不同特征之间的相关性，然后再使用最小二乘法。第二种方法是使用梯度下降法。

3.1.2 梯度下降法

如果把线性回归看作是一个优化问题，那么我们要优化的目标就是损失函数。损失函数是用来衡量样本误差的函数，我们的优化目标是要求得在误差最小的情况下模型参数的值。这里强调一下损失函数和代价函数的区别：
Loss Function(损失函数): the error for single training example;
Cost Function(代价函数)： the average of the loss functions of the entire training set;
线性回归常用的损失函数是均方误差，表达式为：

$$l^{(i)}(\mathbf{w},b)=\frac{1}{2}{\left({\hat{y}}^{(i)}-y^{(i)}\right)}^2$$

$$L(\mathbf{w},b)=\frac{1}{n}\sum _{i=1}^n{l}^{(i)}(\mathbf{w},b)=\frac{1}{n}\sum _{i=1}^n\frac{1}{2}{\left({\mathbf{w}}^{\top }{\mathbf{x}}^{(i)}+b-y^{(i)}\right)}^2$$   其中 $\hat{y}$ 为预测值，$y$ 为真实值
可以直接用公式表达出来的这类解叫作解析解，大多数深度学习模型并没有解析解，只能通过优化算法有限次迭代模型参数来尽可能降低损失函数的值。这类解叫作数值解。

小批量随机梯度下降：先选取一组模型参数的初始值，如随机选取；接下来对参数进行多次迭代，使每次迭代都可能降低损失函数的值。在每次迭代中，先随机均匀采样一个由固定数目训练数据样本所组成的小批量（mini-batch），然后求小批量中数据样本的平均损失和有关模型参数的导数（梯度），最后用此结果与预先设定的学习率的乘积作为模型参数在本次迭代的减小量。
如下式所示：$$(\mathbf{w},b)\leftarrow (\mathbf{w},b)-\frac{\eta }{|\mathcal{B}|}\sum _{i\in \mathcal{B}}{\partial }_{(\mathbf{w},b)}{l}^{(i)}(\mathbf{w},b)$$

学习率($\eta$): 代表在每次优化中，能够学习的步长的大小
批量大小($B$): 是小批量计算中的批量大小batch size
具体内容可以参考https://www.cnblogs.com/lliuye/p/9451903.html

3.2 基于异常检测的线性回归

异常检测中并不会对任何变量给与特殊对待，异常值的定义是基于基础数据点的整体分布，因此需要采用一种更一般的回归建模：即以相似的方式对待所有变量，通过最小化数据对该平面的投影误差确定最佳回归平面。在这种情况下，假设我们有一组变量 $X_1\dots X_d$， 对应的回归平面如下：

$$a_1\cdot X_1+\dots +a_d\cdot X_d+a_{d+1}=0$$

为了后续计算的方便，对参数进行如下约束： $$\sum _{i=1}^d{a}_i^2=1$$   以$L_2$范数作为目标函数： $$L=\left||U\cdot A\right|{|}_2$$

4、主成分分析

4.1 原理

对于 $d$ 维，包含 $N$ 个样本的数据，用 $R_i$ 表示其中第 $i$ 行为：$[x_{i1}...x_{id}]$。由此可以得到 $d\times d$ 的协方差矩阵（标准的PCA应当计算相关系数矩阵，即对数据进行均值为0方差为1的标准化处理，而协方差矩阵只需要减去均值即可）:

$$\Sigma =(R-\bar{R}{)}^T\cdot (R-\bar{R})$$

易知协方差矩阵 $\Sigma$ 是对称并且半正定的，因此可以进行相似对角化：

$$\Sigma =P\cdot D\cdot P^T$$

这里的 $D$ 为对角矩阵，对角元素为特征值；$P$ 为标准正交矩阵，每一行为对应的特征向量；这些标准正交向量提供了数据应该投影的轴线方向。与异常检测相关的主成分分析的主要性质如下：

如果前 $k$ 的特征向量选定之后（根据最大的$k$个特征值），由这些特征向量定义的 $k$ 维超平面是在所有维度为 $k$ 的超平面中，所有数据点到它的均方距离尽可能小的平面。

如果将数据转换为与正交特征向量对应的轴系，则转换后的数据沿每个特征向量维的方差等于相应的特征值。在这种新表示中，转换后的数据的协方差为0。

由于沿特征值小的特征向量的转换数据的方差很低，因此沿这些方向的变换数据与平均值的显着偏差可能表示离群值。

需要注意的是，相比2.2节的内容，这里提供了一个更加普遍的解决方法。2.2中的内容可以归为主成分分析中只保留最大特征值对应的特征向量的情况。

在得到这些特征值和特征向量之后，可以将数据转换到新的坐标系中。以 $Y_1...Y_N$ 表示新坐标系中的数据，这些数据可以通过原始向量 $R_i$ 与包含新轴系的标准正交特征向量矩阵 $P$ 的乘积来实现。 $$Y_i=R_i\cdot P$$

在许多涉及高维数据集的真实场景中，很大一部分特征值往往非常接近于零。这意味着大多数数据都沿着一个低维的子空间排列。从异常检测的角度来看，这是非常方便的，因为离这些投影方向非常远的观测值可以被假定为离群值。例如，对于特征值较小（方差较小）的特征向量 $j$，第 $i$ 条记录的 $y_{ij}$ 与 $y_{kj}$ 的其他值的偏差较大，说明有离群行为。这是因为当$j$固定而$k$变化时，$y_{kj}$ 的值应当变化不大。因此，$y_{ij}$ 值是不常见的。

在不选取任何特定的 $k$ 维集合的情况下，一种更精确的异常检测建模方法是使用特征值来计算数据点沿每个主分量方向到质心的归一化距离。设 $e_j$为第 $j$ 个特征向量，${\lambda }_j$ 为沿该方向的方差(特征值)。数据点$\bar{X}$相对于对数据质心$\bar{\mu} $的总体归一化异常得分可以由下式给出:

$$S\mathrm{core}(\bar{X})=\sum _{j=1}^d\frac{{\left|(\bar{X}-\bar{\mu })\cdot \bar{e}j\right|}^2}{\lambda j}$$

4.2 归一化

两种常用的归一化方法：
一、min-max标准化（Min-Max Normalization）

也称为离差标准化，是对原始数据的线性变换，使结果值映射到[0 - 1]之间。转换函数如下：

其中max为样本数据的最大值，min为样本数据的最小值。这种方法有个缺陷就是当有新数据加入时，可能导致max和min的变化，需要重新定义。
二、Z-score标准化方法

这种方法给予原始数据的均值（mean）和标准差（standard deviation）进行数据的标准化。经过处理的数据符合标准正态分布，即均值为0，标准差为1，转化函数为：

5、作业

使用pyod库生成example并使用该库的pca模块进行检测

import pyod
import pandas as pd 
from pyod.models.pca  import PCA
#pyod.models.hbos.HBOS(alpha=0.1, contamination=0.1, n_bins=10, tol=0.5)
import sklearn
import pyod
import numpy as np
# 训练一个HBOS检测器
clf_name = 'PCA'
clf = PCA() # 初始化检测器clf

# 通过pyod.utils.data.generate_data创建一个toy数据
data = pyod.utils.data.generate_data(n_train=100,
    n_test=50,
    n_features=2,
    contamination=0.1,
    train_only=False,
    offset=10,
    behaviour='old',
    random_state=None,)
X1 = data[0]
X2 = data[2]
X_train = np.r_[X1,X2]



clf.fit(X_train) # 使用X_train训练检测器clf

# 返回训练数据X_train上的异常标签和异常分值
y_train_pred = clf.labels_  # 返回训练数据上的分类标签 (0: 正常值, 1: 异常值)
y_train_scores = clf.decision_scores_  # 返回训练数据上的异常值 (分值越大越异常)

y_train_pred
#array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
   #    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
   #    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
   #    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 #      0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 #      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  #     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1])