内容概要：

1. 引言

真实数据集中不同维度的数据通常具有高度的相关性，这是因为不同的属性往往是由相同的基础过程以密切相关的方式产生的。在古典统计学中，这被称为——回归建模，一种参数化的相关性分析。

回归分析是对具有因果关系的影响因素（自变量）和预测对象（因变量）所进行的数理统计分析处理。只有当自变量与因变量确实存在某种关系时，建立的回归方程才有意义。因此，作为自变量的因素与作为因变量的预测对象是否有关，相关程度如何，以及判断这种相关程度的把握性多大，就成为进行回归分析必须要解决的问题。进行相关分析，一般要求出相关关系，以相关系数的大小来判断自变量和因变量的相关的程度。
——百度百科-回归分析

变量的相关性分析主要分为2类：一类相关性分析试图通过其他变量预测单独的属性值，另一类方法用一些潜在变量来代表整个数据。前者的代表是 线性回归，后者一个典型的例子是 主成分分析。本文将会用这两种典型的线性相关分析方法进行异常检测。

需要明确的是，这里有两个重要的假设：
  假设一：近似线性相关假设。线性相关假设是使用线性回归和主成分分析这两种模型进行异常检测的重要理论基础。

  假设二：子空间假设。子空间假设认为数据是镶嵌在低维子空间中的，线性方法的目的是找到合适的低维子空间使得异常点(o)在其中区别于正常点(n)。

基于这两点假设，在异常检测的第一阶段，为了确定特定的模型是否适合特定的数据集，对数据进行探索性（Exploratory Data Analysis，EDA）和可视化分析是非常关键的。

2. EDA 与 数据可视化

在统计学中，探索性数据分析（EDA）是一种分析数据集以概括其主要特征的方法，通常使用可视化方法。统计模型可用可不用，但主要的是，EDA是为了观察到数据在正式建模或假设检验任务之外的内容。探索性数据分析是John Tukey提出的，鼓励统计学家去研究数据，并尽可能地提出假设，以引领新的数据收集和实验。EDA不同于初始数据分析（Initial Data Analysis，IDA），IDA更集中于检查模型拟合和假设检验所需的假设，并根据需要处理缺失值，以及进行变量转换。EDA包含IDA。
—— 翻译自维基百科

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

数据集来源：

• https://www.ebi.ac.uk/ega/studies/EGAS00001000646
• https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/OPQMVF

基因表达数据反映的是直接或间接测量得到的基因转录产物 mRNA在细胞中的丰度，这些数据可以用于分析哪些基因的表达发生了改变，基因之间有何相关性，在不同条件下基因的活动是如何受影响的。它们在医学临床诊断、药物疗效判断、揭示疾病发生机制等方面有重要的应用。高通量检测基因组 mRNA 丰度的方法主要是 cDNA 微阵列、寡核苷酸芯片，随着 cDNA 微阵列和寡核苷酸芯片等高通量检测技术的发展，我们可以从全基因组水平定量或定性检测基因转录产物 mRNA。

基因表达数据分析的对象是在不同条件下，全部或部分基因的表达数据所构成的数据矩阵。通过对该数据矩阵的分析，可以回答一些生物学问题，例如，基因的功能是什么？在不同条件或不同细胞类型中，哪些基因的表达存在差异？在特定的条件下，哪些基因的表达发生了显著改变，这些基因受到哪些基因的调节，或者控制哪些基因的表达？哪些基因的表达是细胞状态特异性的，根据它们的行为可以判断细胞的状态（生存、增殖、分化、凋亡、癌变或应激等）等等。对这些问题的回答，结合其它生物学知识和数据有助于阐明基因的表达调控路径和调控网络。揭示基因调控路径和网络是生物学和生物信息学共同关注的目标，是系统生物学(SystemsBiology)研究的核心内容。

目前，对基因表达数据的分析主要是在三个层次上进行：1、分析单个基因的表达水平，根据在不同实验条件下，基因表达水平的变化，来判断它的功能，例如，可以根据表达差异的显著性来确定肿瘤分型相关的特异基因。采用的分析方法有统计学中的假设检验等。2、考虑基因组合，将基因分组，研究基因的共同功能、相互作用以及协同调控等。多采用聚类分析等方法。3、尝试推断潜在的基因调控网络，从机理上解释观察到的基因表达数据。多采用反向工程的方法。

—— 基因表达数据的分析 ：https://wenku.baidu.com/view/3089fd3e580216fc700afd08.html

# 导入warnings包，利用过滤器来实现忽略警告语句。
import warnings
warnings.filterwarnings('ignore')
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import missingno as msno

1) 载入训练集和测试集

path = 'dataverse_files/'
f= path+'breast-cancer-unsupervised-ad.csv'
Train_data = pd.read_csv(f)

2) 简略观察数据(head()+shape)

Train_data.head()

5 rows × 31 columns

这份数据没有表头，需要手动添加一下表头。上网查了一下，大致意思是，前面30个特征为通过直接或间接方法测量得到的基因转录产物 mRNA在细胞中的丰度，最后一个特征为是否患病的标签。在标签列中 ，“n”表示“正常”，“o”表示“患病”。

columns = ['gen'+str(i) for i in range(30)]
columns.append('label')
Train_data = pd.read_csv(path, names=columns)
Train_data.head()

5 rows × 31 columns

Train_data.tail()

5 rows × 31 columns

3) 通过describe()来熟悉数据的相关统计量

Train_data.describe()

8 rows × 30 columns

4) 通过info()来熟悉数据类型

Train_data.info()

RangeIndex: 367 entries, 0 to 366 Data columns (total 31 columns):

Column Non-Null Count Dtype

---

0 gen0 367 non-null float64 1 gen1 367 non-null float64 2 gen2 367 non-null float64 3 gen3 367 non-null float64 4 gen4 367 non-null float64 5 gen5 367 non-null float64 6 gen6 367 non-null float64 7 gen7 367 non-null float64 8 gen8 367 non-null float64 9 gen9 367 non-null float64 10 gen10 367 non-null float64 11 gen11 367 non-null float64 12 gen12 367 non-null float64 13 gen13 367 non-null float64 14 gen14 367 non-null float64 15 gen15 367 non-null float64 16 gen16 367 non-null float64 17 gen17 367 non-null float64 18 gen18 367 non-null float64 19 gen19 367 non-null float64 20 gen20 367 non-null float64 21 gen21 367 non-null float64 22 gen22 367 non-null float64 23 gen23 367 non-null float64 24 gen24 367 non-null float64 25 gen25 367 non-null float64 26 gen26 367 non-null float64 27 gen27 367 non-null float64 28 gen28 367 non-null float64 29 gen29 367 non-null float64 30 label 367 non-null object dtypes: float64(30), object(1) memory usage: 89.0+ KB

5)对数值型特征进行相关性分析

numeric_features = columns[:-1]
numeric = Train_data[numeric_features]
correlation = numeric.corr()

f , ax = plt.subplots(figsize = (14, 14))

sns.heatmap(correlation,square = True)
plt.title('Correlation of Numeric Features with Price',y=1,size=16)
plt.show()

6）每个数字特征的分布可视化

f = pd.melt(Train_data, value_vars=numeric_features)
g = sns.FacetGrid(f, col="variable",  col_wrap=6, sharex=False, sharey=False)
g = g.map(sns.distplot, "value", hist=False, rug=True)
plt.savefig('variable.jpg', dpi=300)  # 保存到本地。当图片比较大时，保存到本地再打开会比较清晰。

7）变量两两之间的相关性

sns.set()
# 因为30个特征生成两两之间的相关性图有30x30 共900个子图，子图看起来很密，也比较吃配置。
# 所以这里只展示前10个特征两两之间的相关性。
sns.pairplot(Train_data[numeric_features[:10]],size = 2 ,kind ='scatter',diag_kind='kde')
plt.savefig('correlation.png')
plt.show()

8）数据降维可视化

from sklearn.manifold import TSNE 


tsne = TSNE(n_components=2,  # 嵌入空间的维度（嵌入空间的意思就是结果空间）
            init='pca',  # 嵌入空间的初始化，接收字符串‘random’，‘pca’或者一个numpy数组
            random_state=0)
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]),
                color=plt.cm.Set1(color[label[i]] / 10.),
                fontdict={'weight': 'bold', 'size': 9})
plt.xticks([])
plt.yticks([])
plt.title('Visualization of data dimension reduction')
plt.show()

拓展资料：
TSNE的全称：t-distributed Stochastic Neighbor Embedding
t-SNE（TSNE）将数据点之间的相似度转换为概率。原始空间中的相似度由高斯联合概率表示，嵌入空间的相似度由“学生t分布”表示。
虽然Isomap，LLE和variants等数据降维和可视化方法，更适合展开单个连续的低维的manifold。但如果要准确的可视化样本间的相似度关系，如对于下图所示的S曲线（不同颜色的图像表示不同类别的数据），t-SNE表现更好。因为t-SNE主要是关注数据的局部结构。

• sklearn中tsne可视化

https://www.deeplearn.me/2137.html

• t-SNE 原理及Python实例：

https://zhuanlan.zhihu.com/p/148170862

3. 线性回归

在线性回归中，我们假设不同维度的变量具有一定的相关性，并可以通过一个相关系数矩阵进行衡量。因此对于特定的观测值，可以通过线性方程组来建模。在实际应用中，观测值的数量往往远大于数据的维度，导致线性方程组是一个超定方程，不能直接求解。因此需要通过优化的方法，最小化模型预测值与真实数据点的误差。

线性回归是统计学中一个重要的应用，这个重要的应用往往是指通过一系列自变量去预测一个特殊因变量的值。在这种情况下，异常值是根据其他自变量对因变量的影响来定义的，而自变量之间相互关系中的异常则不那么重要。这里的异常点检测主要用于数据降噪，避免异常点的出现对模型性能的影响，因而这里关注的兴趣点主要是正常值(n)。

而我们通常所说的异常检测中并不会对任何变量给与特殊对待，异常值的定义是基于基础数据点的整体分布，这里我们关注的兴趣点主要是异常值(o)。

**
广义的回归建模只是一种工具，这种工具既可以用来进行数据降噪也可以进行异常点检测。

利用线性回归模型检测异常的关键思想：
真实值与模型预测值之间的差值，用来衡量这个数据点的异常程度。

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

3.1.1 拟合模型的优化——最小二乘法

最小二乘法是由勒让德在19世纪发现的，是一种在误差估计、不确定度、系统辨识及预测、预报等数据处理诸多学科领域得到广泛应用的数学工具。
形式如下式：

观测值就是我们的多组样本，理论值就是我们的假设拟合函数。目标函数也就是在机器学习中常说的损失函数，我们的目标是得到使目标函数最小化时候的拟合函数的模型。举一个最简单的线性回归的简单例子，比如我们有 个只有一个特征的样本：

样本采用一般的为 次的多项式拟合， 拟合函数如下：

其中 为参数。

最小二乘法就是要找到一组 使得预测值和真实值之间的残差平方和最小，

即求 。

这种求解线性回归参数的方法叫最小二乘法（generalized least squares，简称GLS），又称最小平方法（least square method）。最小二乘法是残差满足正态分布的最大似然估计。

拓展：如何由最大似然估计推得最小二乘法：

MLE (最大似然) 与 LS (最小二乘) 与 MAP (最大后验)： https://blog.csdn.net/cdd2xd/article/details/71248570

3.1.2 最小二乘法的直接求解——标准方程法

标准方程法（Normal Equation，又叫正规方程法）是一种基础的最小二乘方法，通过把目标函数转换为矩阵的形式，对矩阵求导来实现最小二乘法的直接求解。

为了简单起见，这里我们一元线性回归为例：

变量Y为因变量，也就是我们要预测的值；为一系列因变量，也就是输入值。系数为要学习的参数。假设数据共包含个样本，第个样本包含的数据为和，带入式(1)如下式所示：

这里为第个样本的误差。以 代表 的因变量矩阵%7D%5E%7BT%7D#card=math&code=%7B%28y%7B1%7D…y%7BN%7D%29%7D%5E%7BT%7D)，即样本中的真实值；以代表#card=math&code=N%20%5Ctimes%20%28d%2B1%29)的自变量矩阵，其中第行为#card=math&code=%28x%7Bj1%7D…x%7Bjd%7D%2C%201%29)；以 代表 的系数矩阵%5E%7BT%7D#card=math&code=%28a%7B1%7D…a%7Bd%2B1%7D%29%5E%7BT%7D)。则模型可表示为：

%20%3D%20U%20%5Ccdot%20A%0A#card=math&code=f%28U%2C%20A%29%20%3D%20U%20%5Ccdot%20A%0A)

定义目标函数为：

%20%3D%20%5Cfrac%7B1%7D%7B2%7D%7B%5Cleft%5C%7C%20%7BY%20-%20U%20%5Ccdot%20A%7D%20%5Cright%5C%7C%5E2%7D%20%0A#card=math&code=L%28A%29%20%3D%20%5Cfrac%7B1%7D%7B2%7D%7B%5Cleft%5C%7C%20%7BY%20-%20U%20%5Ccdot%20A%7D%20%5Cright%5C%7C%5E2%7D%20%0A)

目标函数是关于的凸函数，其对求偏导为：

%7D%7D%7B%7B%5Cpartial%20A%7D%7D%20%3D%20%5Cfrac%7B1%7D%7B2%7D%5Cfrac%7B%7B%5Cpartial%20%7B%7B%5Cleft%5C%7C%20%7BY%20-%20U%20%5Ccdot%20A%7D%20%5Cright%5C%7C%7D%5E2%7D%7D%7D%7B%7B%5Cpartial%20A%7D%7D%20%3D%20%20-%20%7BU%5ET%7D(Y%20-%20U%20%5Ccdot%20A)%0A#card=math&code=%5Cfrac%7B%7B%5Cpartial%20L%28A%29%7D%7D%7B%7B%5Cpartial%20A%7D%7D%20%3D%20%5Cfrac%7B1%7D%7B2%7D%5Cfrac%7B%7B%5Cpartial%20%7B%7B%5Cleft%5C%7C%20%7BY%20-%20U%20%5Ccdot%20A%7D%20%5Cright%5C%7C%7D%5E2%7D%7D%7D%7B%7B%5Cpartial%20A%7D%7D%20%3D%20%20-%20%7BU%5ET%7D%28Y%20-%20U%20%5Ccdot%20A%29%0A)

令%7D%7D%7B%7B%5Cpartial%20A%7D%7D%3D0#card=math&code=%5Cfrac%7B%7B%5Cpartial%20L%28A%29%7D%7D%7B%7B%5Cpartial%20A%7D%7D%3D0)，得到最优参数为：
%5E%7B-1%7D%20%5Ccdot%5Cleft(U%5E%7BT%7D%20%5Ccdot%20Y%5Cright)%0A#card=math&code=A%3D%5Cleft%28U%5E%7BT%7D%20%5Ccdot%20U%5Cright%29%5E%7B-1%7D%20%5Ccdot%5Cleft%28U%5E%7BT%7D%20%5Ccdot%20Y%5Cright%29%0A)

标准方程法要求矩阵 可逆，即是满秩的。当不可逆时可以通过两种方法进行参数估计，一种先使用主成分分析等方法来预处理数据，消除不同特征之间的相关性，然后再使用最小二乘法。第二种方法是使用梯度下降法。

例子：用标准方程法求解参数

# 所求的参数在机器学习回归模型中也叫权重值
def calculate_weights(xArr, yArr):
    xMat = np.mat(xArr)
    yMat = np.mat(yArr)
    xTx = xMat.T * xMat
    # linale.det()矩阵行列式，如果有列或者行全为0，计算结果就为0
    if np.linalg.det(xTx) == 0.0:
        print('This matrix can not do inverse')
        return
    # xTx.I 为 xTx的逆矩阵
    ws = xTx.I * xMat.T * yMat
    return ws


# 使用前面可视化中降维后的数据，这里我们使用标签为“正常”的数据建立模型
sample = result[Train_data['label']=='n']
x_data = sample[:, 0, np.newaxis]
y_data = sample[:, 1, np.newaxis]

# 添加偏置值，代表的是回归方程里面的a_{d+1}
X_data = np.concatenate((np.ones((len(x_data), 1)), x_data), axis=1)

ws = calculate_weights(X_data, y_data)

print(ws)

[[ 0.69108814] [-0.6327987 ]]

使用求解所得的参数代入一元线性回归函数方程里，画出函数的图像。

x_test = [[x_data.min()], [x_data.max()]]
y_test = ws[0] + x_test*ws[1] 

fig = plt.figure(figsize = (7, 7))
plt.scatter(x_data, y_data, color='c', marker='*')
plt.plot(x_test, y_test, color='r')
plt.xticks([])
plt.yticks([])
plt.title('Fitting effect of linear regression', fontsize=15)
plt.show()

参考教程：https://www.bilibili.com/video/BV1Rt411q7WJ?p=15

3.1.2 最小二乘法的辗转求解——梯度下降法

梯度下降法（gradient descent），又名最速下降法（steepest descent）是求解无约束最优化问题最常用的方法，它是迭代法的一种，每一步主要的操作是求解目标函数的梯度向量，将当前位置的负梯度方向作为搜索方向（因为在该方向上目标函数下降最快，这也是最速下降法名称的由来）。

迭代法也称辗转法，是一种不断用变量的旧值递推新值的过程，跟迭代法相对应的是直接法(或者称为一次解法)，即一次性解决问题。前面讲到的标准方程法就是与梯度下降法相对应的直接法。

因此，梯度下降法可以用于求解最小二乘问题（线性和非线性都可以）。

①数据集
监督学习一般靠数据驱动。我们通常收集一系列的真实数据，例如多栋房屋的真实售出价格和它们对应的面积和房龄。我们希望在这个数据上面寻找模型参数来使模型的预测价格与真实价格的误差最小。在机器学习术语里，该数据集被称为训练数据集（training data set）或训练集（training set），通常还应该有一个用于防止过拟合的交叉验证集和一个用于评估模型性能的测试集(test set)。一栋房屋被称为一个样本（sample），其真实售出价格叫作标签（label），用来预测标签的两个因素叫作特征（feature）。

②损失函数
如果把线性回归看作是一个优化问题，那么我们要优化的目标就是损失函数。损失函数是用来衡量样本误差的函数，我们的优化目标是要求得在误差最小的情况下模型参数的值。这里强调一下损失函数和代价函数的区别：

注意：
Loss Function(损失函数)：the error for single training example;
Cost Function(代价函数)：the average of the loss functions of the entire training set;

线性回归常用的损失函数是均方误差，表达式为：

%7D(%5Cmathbf%7Bw%7D%2C%20b)%3D%5Cfrac%7B1%7D%7B2%7D%5Cleft(%5Chat%7By%7D%5E%7B(i)%7D-y%5E%7B(i)%7D%5Cright)%5E%7B2%7D%0A#card=math&code=l%5E%7B%28i%29%7D%28%5Cmathbf%7Bw%7D%2C%20b%29%3D%5Cfrac%7B1%7D%7B2%7D%5Cleft%28%5Chat%7By%7D%5E%7B%28i%29%7D-y%5E%7B%28i%29%7D%5Cright%29%5E%7B2%7D%0A)

%3D%5Cfrac%7B1%7D%7Bn%7D%20%5Csum%7Bi%3D1%7D%5E%7Bn%7D%20l%5E%7B(i)%7D(%5Cmathbf%7Bw%7D%2C%20b)%3D%5Cfrac%7B1%7D%7Bn%7D%20%5Csum%7Bi%3D1%7D%5E%7Bn%7D%20%5Cfrac%7B1%7D%7B2%7D%5Cleft(%5Cmathbf%7Bw%7D%5E%7B%5Ctop%7D%20%5Cmathbf%7Bx%7D%5E%7B(i)%7D%2Bb-y%5E%7B(i)%7D%5Cright)%5E%7B2%7D%0A#card=math&code=L%28%5Cmathbf%7Bw%7D%2C%20b%29%3D%5Cfrac%7B1%7D%7Bn%7D%20%5Csum%7Bi%3D1%7D%5E%7Bn%7D%20l%5E%7B%28i%29%7D%28%5Cmathbf%7Bw%7D%2C%20b%29%3D%5Cfrac%7B1%7D%7Bn%7D%20%5Csum%7Bi%3D1%7D%5E%7Bn%7D%20%5Cfrac%7B1%7D%7B2%7D%5Cleft%28%5Cmathbf%7Bw%7D%5E%7B%5Ctop%7D%20%5Cmathbf%7Bx%7D%5E%7B%28i%29%7D%2Bb-y%5E%7B%28i%29%7D%5Cright%29%5E%7B2%7D%0A)

其中 为预测值， 为真实值。

③优化算法 - 随机梯度下降
当模型和损失函数形式较为简单时，上面的误差最小化问题的解可以直接用公式表达出来。这类解叫作解析解（analytical solution）。本节使用的线性回归和平方误差刚好属于这个范畴。然而，大多数深度学习模型并没有解析解，只能通过优化算法有限次迭代模型参数来尽可能降低损失函数的值。这类解叫作数值解（numerical solution）。

在求数值解的优化算法中，小批量随机梯度下降（mini-batch stochastic gradient descent）被广泛使用。它的算法很简单：先选取一组模型参数的初始值，如随机选取；接下来对参数进行多次迭代，使每次迭代都可能降低损失函数的值。在每次迭代中，先随机均匀采样一个由固定数目训练数据样本所组成的小批量（mini-batch），然后求小批量中数据样本的平均损失和有关模型参数的导数（梯度），最后用此结果与预先设定的学习率的乘积作为模型参数在本次迭代的减小量。如下式所示：

%20%5Cleftarrow(%5Cmathbf%7Bw%7D%2C%20b)-%5Cfrac%7B%5Ceta%7D%7B%7C%5Cmathcal%7BB%7D%7C%7D%20%5Csum%7Bi%20%5Cin%20%5Cmathcal%7BB%7D%7D%20%5Cpartial%7B(%5Cmathbf%7Bw%7D%2C%20b)%7D%20l%5E%7B(i)%7D(%5Cmathbf%7Bw%7D%2C%20b)%0A#card=math&code=%28%5Cmathbf%7Bw%7D%2C%20b%29%20%5Cleftarrow%28%5Cmathbf%7Bw%7D%2C%20b%29-%5Cfrac%7B%5Ceta%7D%7B%7C%5Cmathcal%7BB%7D%7C%7D%20%5Csum%7Bi%20%5Cin%20%5Cmathcal%7BB%7D%7D%20%5Cpartial%7B%28%5Cmathbf%7Bw%7D%2C%20b%29%7D%20l%5E%7B%28i%29%7D%28%5Cmathbf%7Bw%7D%2C%20b%29%0A)

学习率(): 代表在每次优化中，能够学习的步长的大小
批量大小(): 是小批量计算中的批量大小batch size

# 计算平均误差
def compute_error(x_data, y_data, b, k):
    """
    x_data:样本特征值
    y_data:样本标签值
    b:拟合回归线的截距
    k:拟合回归线的斜率
    return：平均误差平方和
    """
    x_data, y_data = np.array(x_data), np.array(y_data)
    total_error = np.sum(np.square(y_data - (x_data*k+b)))
    return total_error / len(x_data) / 2


# 梯度下降法
def gradient_descend_runner(x_data, y_data, b, k, lr, epochs):
    """
    Δθ0 = 1/m * Σ(h(xi)-yi) * lr
    Δθ1 = 1/m * Σ[(h(xi)-yi)xi] * lr
    θ0 = θ0 - Δθ0
    θ1 = θ1 - Δθ1
    """
    x_data, y_data = np.asarray(x_data), np.asarray(y_data)
    m = len(x_data)
    for i in range(epochs+1):
        "迭代到一定次数停止"
        b_grad = np.sum((x_data*k+b) - y_data) / m
        k_grad = np.sum(((x_data*k+b) - y_data) * x_data) / m
        b = b - b_grad * lr
        k = k - k_grad * lr
    if i%30000 == 0 :
        print('epochs: ', i)
        print('b = {0}, k = {1}, error = {2}'.format(b, k, compute_error(x_data, y_data, b, k)))
        plt.plot(x_data, y_data, 'b.')
        plt.plot(x_data, x_data*k+b, 'r')
        plt.show()

    return b, k

b = 0
k = 0
lr = 0.001
epochs = 100000  # 迭代周期

print('Starting...b={}, k={}, error={}'.format(b, k, computer_error(x_data, y_data, b, k)))

b, k = gradient_descend_runner(x_data, y_data, b, k, lr, epochs)
print('Finished...b={}, k={}, error={}'.format(b, k, computer_error(x_data, y_data, b, k)))

# 画出回归线
plt.figure(figsize=(6,9))
plt.scatter(x_data, y_data, color='c', marker='*')
plt.plot(x_data, x_data*k+b, color='r')
plt.xlim([0.,0.6])
plt.ylim([0.,0.9])
plt.xticks([])
plt.yticks([])
plt.title('Fitting effect after 100000 epochs using gradient descent',
          fontsize=15)
plt.show()

Starting…b=0, k=0, error=0.1590259669541645

epochs: 0
b = 0.0005493323061646533, k = 0.00011288700050332632, error = 0.15871161239153864

epochs: 30000
b = 0.5830887810709732, k = -0.15806236274970362, error = 0.006737432559999097

epochs: 60000
b = 0.6228294568495886, k = -0.3327517630105567, error = 0.005648881447415392

epochs: 90000
b = 0.647946700101129, k = -0.44316045809927257, error = 0.005214048700840199

Finished…b=0.6540649743193592, k=-0.4700547593023905, error=0.00513782721607625

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

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

为了后续计算的方便，对参数进行如下约束：

以范数作为目标函数：

这样的一个问题可以通过主成分分析方法得到有效解决，我们会单独用一个部分进行讨论。

4. 主成分分析

上一节的最小二乘法试图找到一个与数据具有最佳匹配 维超平面。主成分分析方法可用于解决这一问题的广义版本。具体来说，它可以找到任意 #card=math&code=k%28%20k%3Cd%20%29) 维的最优表示超平面，从而使平方投影误差最小化。

4.1 原理推导

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

%5E%7BT%7D%20%5Ccdot%20(R%20-%20%5Cbar%7BR%7D)%20%0A#card=math&code=%CE%A3%20%3D%20%28R%20-%20%5Cbar%7BR%7D%29%5E%7BT%7D%20%5Ccdot%20%28R%20-%20%5Cbar%7BR%7D%29%20%0A)

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

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

• 如果前 的特征向量选定之后（根据最大的个特征值），由这些特征向量定义的 维超平面是在所有维度为 的超平面中，所有数据点到它的均方距离尽可能小的平面。
• 如果将数据转换为与正交特征向量对应的轴系，则转换后的数据沿每个特征向量维的方差等于相应的特征值。在这种新表示中，转换后的数据的协方差为0。
• 由于沿特征值小的特征向量的转换数据的方差很低，因此沿这些方向的变换数据与平均值的显着偏差可能表示离群值。

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

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

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

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

%3D%5Csum%7Bj%3D1%7D%5E%7Bd%7D%20%5Cfrac%7B%5Cleft%7C(%5Cbar%7BX%7D-%5Cbar%7B%5Cmu%7D)%20%5Ccdot%20%5Cbar%7Be%7D%7Bj%7D%5Cright%7C%5E%7B2%7D%7D%7B%5Clambda%7Bj%7D%7D%0A#card=math&code=S%20%5Coperatorname%7Bcore%7D%28%5Cbar%7BX%7D%29%3D%5Csum%7Bj%3D1%7D%5E%7Bd%7D%20%5Cfrac%7B%5Cleft%7C%28%5Cbar%7BX%7D-%5Cbar%7B%5Cmu%7D%29%20%5Ccdot%20%5Cbar%7Be%7D%7Bj%7D%5Cright%7C%5E%7B2%7D%7D%7B%5Clambda%7Bj%7D%7D%0A)

值得注意的是，对异常得分的大部分贡献是由 值较小的主成分的偏差提供的，这一点上文中有提及过。主成分分析比因变量回归能更稳定地处理少数异常值的存在。这是因为主成分分析是根据最优超平面来计算误差的，而不是一个特定的变量。当数据中加入更多的离群点时，最优超平面的变化通常不会大到影响离群点的选择。因此，这种方法更有可能选择正确的异常值，因为回归模型一开始就更准确。

4.2 归一化问题

当不同维度的尺度差别较大时，使用PCA有时并不能得到直观有效的结果。例如，考虑一个包含年龄和工资等属性的人口统计数据集。工资属性的范围可能是几万，而年龄属性几乎总是小于100，使用主成分分析会导致主成分被高方差属性所控制。对于一个只包含年龄和工资的二维数据集，最大的特征向量几乎与工资轴平行，这会降低异常点检测过程的有效性。因此，一个自然的解决方案是对数据进行均值为0方差为1的标准化处理。这隐含地导致在主成分分析中使用相关矩阵而不是协方差矩阵。当然，这个问题并不是线性建模所独有的，对于大多数异常检测算法，都需要使用这样的预处理。

5. 回归分析的局限性

回归分析作为检测离群值的工具有一些局限性。这些缺点中最重要的是在本章的一开始就讨论了，其中探讨了回归分析的数据特定性质。特别是，为了使回归分析技术有效，数据需要高度相关，并沿着低维子空间对齐。当数据不相关，但在某些区域高度聚集时，这种方法可能不会有效。

另一个相关的问题是，数据中的相关性在本质上可能不是全局性的。最近的一些分析观察表明，子空间相关性是特定于数据的特定位置的。在这种情况下，由主成分分析发现的全局子空间对于异常检测是次优的。因此，为了创建更一般的局部子空间模型，有时将线性模型与邻近模型(在后续章节中讨论)结合起来是有用的。这将是高维和子空间异常检测的主题，将在后续章节详细讨论。

6. 总结

真实数据中，数据不同属性之间往往具有显著的相关性。在这种情况下，线性建模可以提供一种有效的工具来从底层数据中移除异常值或者进行异常检测。对于其他基于因变量回归的应用，线性建模是一种工具，去除异常值对于提高此类应用的性能是非常重要的。在大多数情况下，主成分分析提供了去除异常值和进行异常检测最有效的方法，因为它对存在少数异常值的数据更有鲁棒性。

7. 练习

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

contamination = 0.1  # 异常值的比例
n_train = 20000  # 训练集样本容量
n_test = 10000  # 测试样本容量

# 生成样本数据集
X_train, y_train, X_test, y_test =  generate_data(n_train=n_train, n_test=n_test,
                                                   n_features=200,  contamination=contamination,
                                                   random_state=2021)
# 模型训练
model_name = 'PCA'
model = PCA(n_components=2)  # 子空间的维度，这里我们把20维的数据降到2维。
model.fit(X_train)

y_train_pred = model.labels_  # 训练集训练后生成的标签，0为正常值，1为异常值
y_train_scores = model.decision_scores_   # 训练数据的异常值得分
params = model.get_params()  # 模型参数的估计量

# 对测试集进行预测
y_test_pred = model.predict(X_test)  # 生成测试集的预测标签
y_test_scores = model.decision_function(X_test)  # 使用合适的探测器预测X的原始异常分数。

# 打印评估结果
# print("\nParams: ", params)
print("\nOn Training Data:")
evaluate_print(model_name, y_train, y_train_scores)
print("\nOn Test Data:")
evaluate_print(model_name, y_test, y_test_scores)

返回结果：

Params: {‘contamination’: 0.1, ‘copy’: True, ‘iterated_power’: ‘auto’, ‘n_components’: 2, ‘n_selected_components’: None, ‘random_state’: None, ‘standardization’: True, ‘svd_solver’: ‘auto’, ‘tol’: 0.0, ‘weighted’: True, ‘whiten’: False}
On Training Data: PCA ROC:1.0, precision @ rank n:1.0
On Test Data:

PCA ROC:1.0, precision @ rank n:1.0

参考资料：

• pyod.models.pca API文档：

https://pyod.readthedocs.io/en/latest/pyod.models.html#module-pyod.models.pca

• 异常检测学习资源（Anomaly Detection Learning Resources）:

https://gitee.com/zhangjxcasm/anomaly-detection-resources

参考文献

[1]:Outlier Analysis Charu C. Aggarwal

[2]:Anomaly Detection: A Survey VARUN CHANDOLA, ARINDAM BANERJEE, and VIPIN KUMAR University of Minnesota

[3]:Outlier Analysis Charu C. Aggarwal

[4]:Anomaly Detection: A Tutorial

[5]:Data Mining Concepts and Techniques Third Edition
**