如何说明回归模型效果好 回归模型的优缺点

机器学习

• 五、基础机器学习方法

• 1.线性回归（Linear Regression）

• 1.1 线性回归的定义
• 1.2 最大似然估计（Maximum Likelihood Estimation）
• 1.3 线性回归的过拟合
• 1.4 梯度下降
• 1.5 梯度下降在python中的应用

• 2.主成分分析（Principal Component Analysis）

• 2.1 主成分分析的数学求解
• 2.2 主成分分析中的投射（Objective）
• 2.3 高维空间下的PCA
• 2.4 主成分分析的主要步骤
• 2.5 主成分分析的python应用

• 3. 聚类（Clustering）

• 3.1 聚类分析的简介
• 3.2 聚类中的距离和相似度
• 3.3 分层聚类
• 3.4 K均值算法

• 4.高斯混合模型（Gaussian Mixture Models）

• 4.1 高斯混合模型（GMM）的由来
• 4.2 最大似然估计求解GMM
• 4.3 EM算法
• 4.4 高斯混合模型的优缺点
• 4.5 高斯混合模型的python应用

五、基础机器学习方法

1.线性回归（Linear Regression）

回归问题要解决的是通过观测数据及其噪音，推断出一个可以描述数据变化趋势的模型，这个模型可以用来估计并预测未知数据。在回归问题中，输出可取任意连续的实数。如下图所示，回归模型包含了已经生成了数据（蓝色曲线），并在相应的输入处指示了功能值的测量噪声（橙色分布）。

回归问题属于有监督学习，回归模型其实是数据点和对应的目标值之间的函数关系。而线性回归模型就是当模型的函数为线性时的回归问题。

1.1 线性回归的定义

线性回归（linear regression）的似然函数（likelihood function）定义为下式：

它的函数y = x^Tθ + ε呈现一条直线，似然概率p(y|x, θ) = N(y|x^Tθ, σ²) 是在x^Tθ处评估的y的概率密度函数。唯一的不确定性源于观察噪声ε（假设ε是已知的）。

给出训练集{(x₁,y₁), … , (x_N,y_N)}，其中y是目标值（标签）。每个数据直接是有条件的线性独立（conditionally independent ）。利用这个训练集可以估计出线性回归模型的参数（Parameter estimation）。

1.2 最大似然估计（Maximum Likelihood Estimation）

模型表述的是每个点位置可能出现的似然概率，当目标值正好在这个似然概率最大最时取到，说明这个模型可以最真实的反应数据和目标值之间的关系。所以可能通过使得似然概率最大来找到模型的最佳参数。

上述n个数据点的总的似然概率是每个点的似然概率乘积，即

最大化似然意味着在给定模型参数的情况下最大化训练数据的预测分布。 我们获得最大似然参数为

对似然概率求log之后再变成负，这样求解最大似然概率变成了最小负log似然概率（negative log-likelihood）

带入线性回归的模型方程y = x^Tθ + ε，而且因为噪音ε的似然概率为高斯分布（假设其标准差为σ），上式中的log p()部分就是（其中const代表常数）

将常数部分分离出并用矩阵X和y代替x_N和y_N得到

上述式子就是从概率角度找到的最优化的损失函数。从统计的角度，损失函数是目标值和预测值的差的平方。这两个损失函数式子是等价的。从这两个式子可以求出模型的参数

实际中的数据分布很难呈现直线性，一条直线很难拟合成复杂的数据集。这时要把线性回归的定义扩大，线性是体现在参数上（“linear in the parameters”），而不是整体一条直线。即把单一的数字x换成一个函数φ(x)。这样模型所代表的是非线性的关系

我们将其中的φ称为k维的特征向量（feature vector），它是输入x的非线性转化。而模型的参数θ依旧是线性的。所有N个φ组成一个特征矩阵Φ（feature matrix），矩阵的维度是N*k。之前的损失函数中的x输入矩阵也可以转化为Φ,损失函数变成

对其求解模型参数θ得到

1.3 线性回归的过拟合

当线性回归方程中的参数过于复杂或者多元线性回归中的元数过多，都可能出现过拟合现象，即低的训练集误差但是有高的测试集误差。为了减小过拟合，在损失函数后面加上一项正则项，使得整个参数矩阵变得稀疏（sparser），减小模型的复杂度。加入L1正则化的回归模型称为Lasso回归（Lasso Regression），加入L2正则化的回归模型称为Ridge回归（Ridge Regression）。这两个模型的具体应用见机器学习从入门到死亡(上) 一文中的第四章1.8节。

1.4 梯度下降

使用1.2节中的方程求解参数θ有两个条件：一是Φ^TΦ是可逆矩阵，二是Φ^TΦ的维度不可以太大，否则求解难度和计算复杂度很大。当输入特征矩阵Φ的维度较大，即数据数量较多时，可以使用梯度下降方法（gradient descent）求解模型参数θ。

使用1.2节中的方程求解参数相当于求解一阶导数为0时的极值点。除了求导的计算难度大，还有可能面临局部最小值甚至最大值问题。因为函数的极值点可能是极大值点。针对这种情况可以通过求解二阶导数，判断极大值或极小值。如果二阶导数小于0，该极值点为极大值。如果二阶导数大于0，该极值点为极小值。

而对于求导的计算难度，使用梯度下降方法代替直接求导。相当于从函数曲线的某一个点开始，按照这个点梯度最大的反方向（负梯度方向），减小一个数值（步长step），到了一个更小的点，继续这么做，直到新的点的梯度为水平，即达到最低点（或者局部最低点），这个点就是梯度下降求到的最小值（或者局部最小值）。如下图所示，展示了从-6点，-2点和-1点，2点分别进行梯度下降，最后达到了两个极小值点。

梯度下降方法是一阶优化算法。梯度下降方法需要减去与该函数在当前点的梯度负值成比例的步长（step），从而找到函数的局部最小值。梯度下降方法利用了：如果一个点从x₀朝向f(x₀)函数在x₀处的负梯度方向即 -(∇f(x₀))^T移动，则f(x₀)下降最快。下降一个步长得到更小的点x₁。使用函数表示为（其中γ为步长step size，要大于0）

重复这一步骤，直到找到局部最小值或者达到收敛或阈值，就可以利用梯度下降方法求解损失函数最小时的模型参数。

梯度下降方法有一个问题就是当数据量大时，求解训练集中所有数据点的梯度用时增多，需要每次下降的求解梯度计算增多，总的时间耗费较大。改进方法是使用随机梯度下降（Stochastic Gradient Descent）。随机梯度下降使用一小部分随机数据（smaller minibatch）代替整体求解梯度，这样每一次下降的计算量减小，用时减小。如果梯度下降要使用训练集中N个点计算梯度导致计算量大，那么随机梯度下降仅从训练集中随机选择M个数据（M远小于N），这样每次的计算就减少了。

梯度下降方法在随机初始化参数θ后使用下面式子更新θ：

而随机梯度下降方法在随机初始化参数θ后使用下面式子更新θ：

梯度下降中除了参数的初始化，还有一个超参数，步长。步长（step-size）也可以叫学习速率（(learning rate），对梯度下降的结果有着重大影响。如果步长太小，每次下降的幅度很小，学习速率低，整个梯度下降的耗时过大。如果步长太大，下降过程中可能直接跳过了极小值点，无法收敛（converge）甚至发生diverge，也就无法找到局部最小值。当函数值在梯度阶跃（下降）之后增加时，说明步长太大。；应该撤消这步并减小步长。当函数值在梯度阶跃之后减小，说明步长可以采用更大值；可以 尝试增加步长。当选择步长时，应该从大到小尝试， 比如第i步的步长γ_i =1/(i + 1)。

如果优化曲面的曲率使得某些区域的缩放比例较差，则梯度下降的收敛可能会非常慢。为了加快收敛速度，将梯度下降改进成带动量的梯度下降（Gradient Descent With Momentum）。具有动量的梯度下降是一种引入附加动量或者记忆量以记住上一次迭代中发生的情况。该记忆量可抑制振荡并平滑梯度更新，从而使得梯度下降保持一个较为平均的速率，减少梯度下降的更新次数，使得收敛加快。下面式子中的∆x_i就是每一次迭代中添加的之前的动量，α∈[-1,1]是动量的系数。

1.5 梯度下降在python中的应用

使用随机梯度下降（SGD）求解线性回归（可以带正则化）的参数，在python的sklearn库中对应sklearn.linear_model.SGDRegressor函数。该函数通过使用随机梯度下降和最小化正则经验损失来拟合线性回归模型。函数的输入要求浮点值化的密集numpy数组的数据。函数具体为sklearn.linear_model.SGDRegressor(loss='squared_loss', *, penalty='l2', alpha=0.0001, l1_ratio=0.15, fit_intercept=True, max_iter=1000, tol=0.001, shuffle=True, verbose=0, epsilon=0.1, random_state=None, learning_rate='invscaling', eta0=0.01, power_t=0.25, early_stopping=False, validation_fraction=0.1, n_iter_no_change=5, warm_start=False, average=False)，其中loss是损失函数的选择，默认为最小二乘；penalty是惩罚项，用于添加l1或l2或elasticnet正则化，默认是l2正则化，elasticnet表示l1和l2的结合；alpha是正则化的系数，alpha值越高，正则化越强；l1_ratio是当正则化选择elasticnet时Elastic Net混合l1正则化的比例，从0到1；fit_intercept代表是否应该估计截距；max_iter是训练数据的最大迭代次数；tol代表停止标准（停止的阈值），如果不是None，则连续n_iter_no_change个时期（损失> best_loss-tol）时训练将停止；early_stopping代表早停止选项，当验证分数没有提高时，是否使用提前停止来终止训练；random_state设置为True时，用于对数据进行随机播放；learning_rate设置学习率；average设置为True时，将计算所有更新的平均SGD权重并将结果存储在coef_属性中，如果将int设置为大于1，则一旦看到的样本总数达到平均值，就会开始平均。具体应用如下：

import numpy as np #导入numpy库用于随机数生成
from sklearn.linear_model import SGDRegressor #导入SGD回归函数
from sklearn.pipeline import make_pipeline #用于数据生成
from sklearn.preprocessing import StandardScaler #用于预处理
n_samples, n_features = 10, 5
rng = np.random.RandomState(0) #用于生成随机数作为训练集
y = rng.randn(n_samples) #生成目标值
X = rng.randn(n_samples, n_features)
#始终缩放输入，最方便的方法是使用pipeline.
reg = make_pipeline(StandardScaler(),SGDRegressor(max_iter=1000, tol=1e-3)) #将数据预处理和SGD回归模型结合生成reg
reg.fit(X, y) #向reg模型输入数据
#结果是 Pipeline(steps=[('standardscaler', StandardScaler()),('sgdregressor', SGDRegressor())])

#调用reg.predict(X)可以使用线性模型进行预测
#调用reg.score(X, y[, sample_weight])返回预测的相关性系数coefficient

随机梯度下降除了用于线性回归问题，还可以用于分类问题。使用随机梯度下降求解出分界线的参数，然后用于分类。SGD的线性分类器在python的sklearn库中对应sklearn.linear_model.SGDClassifier函数。该函数的输入数据表示为要素的浮点值的密集或稀疏数组的数据。默认情况下，该函数适用于线性支持向量机（SVM）。函数具体为sklearn.linear_model.SGDClassifier(loss='hinge', *, penalty='l2', alpha=0.0001, l1_ratio=0.15, fit_intercept=True, max_iter=1000, tol=0.001, shuffle=True, verbose=0, epsilon=0.1, n_jobs=None, random_state=None, learning_rate='optimal', eta0=0.0, power_t=0.5, early_stopping=False, validation_fraction=0.1, n_iter_no_change=5, class_weight=None, warm_start=False, average=False)，其中loss代表的损失函数默认为hinge函数，其余参数与上面的SGDRegressor函数一直。具体应用是：

import numpy as np
from sklearn.linear_model import SGDClassifier
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
X = np.array([[-1, -1], [-2, -1], [1, 1], [2, 1]])
Y = np.array([1, 1, 2, 2]) #分类的标签
clf = make_pipeline(StandardScaler(),SGDClassifier(max_iter=1000, tol=1e-3))
clf.fit(X, Y)
#clf模型的结果是 Pipeline(steps=[('standardscaler', StandardScaler()),('sgdclassifier', SGDClassifier())])
print(clf.predict([[-0.8, -1]]))
#分类结果是[1]

2.主成分分析（Principal Component Analysis）

在处理高纬度数据（比如照片）或者多成分（factors）数据时，为了减少计算量和存储空间，可以对数据进行降维处理，提取其中的关键信息。高纬度数据一般都是高重复数据（overcomplete），数据蕴含的信息有重复信息和多余信息（redundant），其中一些维度上的信息是可以通过其他维度信息的组合得到的。而且高维数据中一些维度上的信息是非重要信息，它对整体数据的分布存在较小的影响。所以从高维数据中提取一些重要维度上的数据，可以在减小计算的同时保证信息损失最小。比如下图，左边的图中信息有x1和x2两个维度；右边的信息就只有x1一个维度。

可是左边的图中在x2的方向上并没有太多信息，而右边图中x1方向的数据依旧反应出整体数据在主要方向上的分布状态。从左边图变化到右边的图就是从二维数据降维到了一维数据，而且最大程度保证了信息的完整。降维环节中最重要的一步就是找到可以反应最多信息的主要方向，使得降维的损失最小化，而这个主要方向就是主成分。找到最重要的几个主成分，并将它们的信息提取出来以此降维的方法就是主成分分析（Principal Component Analysis，简称PCA）。

2.1 主成分分析的数学求解

原始训练集X={x₁, … , x_N}，x_n∈R^D，假设数据经过了预处理，它的均值（mean）为0，数据集的协方差矩阵（covariance matrix）S是

将X乘以一个映射矩阵B压缩到一个低纬度的矩阵Z，其中Z中的每一个数据z_n是

而这个映射矩阵（projection matrix）B=[b₁, … , b_M]∈R^D*M。PCA过程所做的是，先找到原始数据X的压缩版本Z。 然后将压缩后的数据Z重构为X‘，X’的数据位于原始数据空间中，但是具有比X固有的低维表示形式。即从X降维到了X‘。其中重构过程就是压缩矩阵Z乘回映射矩阵B，即X’=B * z=B * B^T * x。映射矩阵B又称为解码器（decoder），B^T称为编码器（encoder）。整个压缩与解压缩过程如下图所示

PCA过程中如何找到一个合适的压缩版本Z，是根据信息损失量来寻找的。PCA要使得数据在低纬度表示时的方差最大化，并以此保留做多的信息。在数据压缩后，保留的大多数信息数量等效于在低维表示中捕获最大的方差（Retaining most information after data compression is equivalent to capturing the largest amount of variance in the low-dimensional code. 源自Hotelling, 1933年提出）。压缩矩阵z的方差是（式子中看出方差和均值无关，所以可以假设z的均值为0）

为了使得低纬度矩阵的方差最大，可以先从映射矩阵中寻找一个向量b₁，使得它映射之后的数据方差先到最大值。这也是寻找压缩矩阵z中的第一个坐标方向z₁，并使得z₁的方差V₁最大。

将上面式子展开并用协方差矩阵S表示为

为了便于计算，可以通过数据平移缩放使得b₁的方差为1，在这个基础上求解上面式子的最大值情况下的b₁。借助数学中的拉格朗日定理（Lagrangian，这里不做过多推导）

并对上式子分别求b₁和λ的偏导数，得到

求解最大值也就是要导数为0，所以上面两个偏导数为0得到

结合特征向量和特征值的定义可以看出b₁是协方差矩阵S的特征向量（eigenvector），λ₁是b₁对应的特征值（eigenvalue）。那么方差V₁求解完就等于λ₁。所以投影到一维子空间上的数据的方差等于与跨越该子空间的特征向量b₁对应的特征值。那么为了最大化低维数据的方差，应该选择与数据的协方差矩阵的最大特征值相关的特征向量，并把该特征向量被称为第一主成分量（the first principal component）。在提取出第一个主成分之后，可以继续提取所需的m个主成分，将原始数据从D维降维到了M维。新的主成分提取和之前的主成分相关，可以使用递归的思想依次提取。假设依旧找到了协方差矩阵S的m-1个特征值和特征向量，现在去找第m个主成分。假设前m-1个主成分分别是b₁到b_m-1，剩余的信息量（remaining

information）可以构成一个新的数据矩阵$\hat{X}$（X是原来的数据矩阵）

$\hat{X}$对应的协方差矩阵是$\hat{S}$，那么第m个主成分也就是最大化$\hat{X}$的方差，如下式：

第m个主成分b_m是$\hat{S}$的最大特征值所对应的特征向量，而且满足下面式子

所以b_m不仅是$\hat{S}$的特征向量，也是S的特征向量。λ_m是$\hat{S}$最大的特征值，它对应也是S的第m大特征值，其在S和$\hat{S}$所对应的特征向量都是b_m。对于之前m-1个主成分，b₁到b_m-1它们也是$\hat{S}$的特征向量，只不过它们的特征值为0。通过上面式子，结合下面式子

可以得到数据映射到第m个主成分的方差是

这意味着，当投影到M维子空间上时，数据的方差等于与数据协方差矩阵的对应特征向量关联的特征值之和。即降维到M维时数据包含的信息即它的方差是（λ_m是协方差矩阵S的第m大特征值）

对于降维到M维时损失的信息（即损失数据的方差）是

所以求解压缩矩阵Z，要先求解出数据X的协方差矩阵S的前m个最大特征值对应的特征向量，将它们作为前m个主成分，然后将数据投射到m个主成分上进行压缩。

2.2 主成分分析中的投射（Objective）

PCA也可以视作直接最小化平均重建误差（the average reconstruction error）的算法。因为投射之后的新数据矩阵$\hat{X}$要和之前的X相似即保留同样信息，所以X和$\hat{X}$的差值（欧式距离）也要最小化。对于数据向量x所在向量空间x∈R^D中，找到一组正交基底（orthonormal basis）(b₁, … , b_D)，x用这组基底表示为

压缩之后的$\hat{x}$是在R^D的一个低纬度子空间中的向量，而且新的维度为M，M小于D，$\hat{x}$用坐标和基底表示为（其中(b₁, … , b_M)也是正交基底，而空间U也叫做主要子空间principal subspace）

为了使得X和$\hat{X}$的差值（欧式距离）也要最小化，可以先把这个差值定义出来，即欧式距离的平均方差（the average squared Euclidean distance），也称为reconstruction error，使用数学表达式J_M为

下面就是寻找合适主要子空间的正交基底和其对应的坐标z_n，使得上式最小。先用导数求解坐标，再寻找最佳基底。对J_M求解z_in的偏导数（partial derivatives）为

根据求导的乘法法则将上面第二个和第三个式子代入第一个式子得到

当导数为0时J_M最小，所以最佳坐标z_in的解为

从上面式子可以看出投影$\hat{x}$的最佳坐标z_in是原始数据点x_n在b_i跨越的一维子空间上的正交投影的坐标。原始数据x的最佳线性投影$\hat{x}$也是正交投影。而在基底(b₁, … , b_M)上的$\hat{x}$的坐标是原始数据x_n在主子空间上的正交投影的坐标。

从线性代数的知识可以得到如果(b₁, … , b_D)是正交基底，那么x在该D维空间的投影$\hat{x}$满足

扩展到这个D维空间的一个M维的子空间上，投影$\hat{x}$满足（其中B是新的基底(b₁, … , b_M)）

虽然$\hat{x}$仍然有D个坐标，但是只需用其中M个坐标表示$\hat{x}$，而其余D-M个坐标映射到基底(b_M+1, … , b_D)上则是0。到此就找到了在主要子空间上最佳投射坐标$\hat{x}$，下面继续求解这个最佳基底（best basis）。从之前$\hat{x}$的坐标可以得到

原始数据x_n可以使用所有D个基底向量b_d的线性结合（linear combination）表示并将其分解成D个主成分向量和D-M个非主成分向量

x_n和它的投影$\hat{x}$_n的差值就是原始数据点和投影的不同，两者做减法得到差值向量（displacement vector）

对这个差值向量求模再求平均（平均平方重构误差）并将它定义为投影之后的损失函数，也就是先前提到的J_M，这里也验证了之前设置的欧式距离作为差值推出的损失函数，两者相等。

对J_M式子拆分后化简并结合trace的定义得到（其中S就是协方差矩阵）

将损失函数公式化为数据的协方差矩阵，并将其投影到主要子空间的正交互补上。因此，最小化平均平方重构误差等效于最小化投影到被忽略的子空间（即主子空间的正交补数）时的数据方差。同样，将保留在主要子空间中的投影的方差最大化，投影损失函数也就化为在上一节中PCA的最大方差公式。即当投影到M维主子空间时，平均平方重构误差为（其中λ_j是协方差矩阵的特征值，这个式子也和上一节中的PCA数据损失方差一样了）

为了使上面损失函数最小化，需要选择最小的D-M特征值。 它们对应的特征向量是主要子空间正交互补的基底。这意味着主要子空间的基底就是数据协方差矩阵的最大M个特征值相关联的特征向量(b₁, … , b_M)。至此，从投射的角度推导出使用PCA求解最佳基底。

2.3 高维空间下的PCA

为了求解PCA，需要先求解协方差矩阵S，S是D * D的矩阵。但是在高维空间，维度D很大，求解S计算量庞大。为了减小计算量，使用一个求解等价特征向量的方法替换直接求解S。对于原始数据点，一般数据点的个数会远小于维度，即N远小于D。S是X * X^T，是D * D维，但是如果可以使用X^T * X进行替代，维度就减小为N * N，计算量大幅减小。而事实上，X * X^T和X^T * X两个矩阵的特征向量是一样的。

PCA中有下面等式

从协方差矩阵S定义可以得到

将上面式子左边和右边都左乘（left-multiply）以X^T得到

从这个式子右边的化简后部分可以看出，X * X^T和X^T * X两个矩阵的非零特征值是一样的，都是λ_m。而1/N * X^T * X的λ_m对应的特征向量c_m是

对比1/N * X * X^T的特征向量，可以得到两个矩阵的非零特征向量也是一样的。此时计算的X^T * X只是一个N * N的矩阵，计算量大幅减小了。当将特征向量等式同时左乘X得到

这个式子的左边括号内部分就是协方差矩阵S，所以可以得到X * c_m也是S的一个特征向量。至此可以完全使用X^T * X代替X * X^T求解PCA过程中的特征向量和特征值。

2.4 主成分分析的主要步骤

主成分分析主要包含五步：第一步将原始数据同时减去均值（Mean subtraction），将数据集的中心移到原点，此时新的数据集的均值为0；第二步将数据标准化（Standardization），对于D个维度的每一个维度将数据除以标准差，此时新的数据集在每一个轴（axis）上的方差为1；第三步是求解协方差矩阵bingo对其进行特征分解（Eigendecomposition），得到特征向量和特征值，大的特征值对于主要子空间（principal subspace）；第四步将标准化之后的原始数据投影（Projection）在主要子空间上，得到投射之后的坐标z和投射之后数据集$\hat{x}$；第五步将$\hat{x}$进行反向去除标准化（undo the standardization），乘以标准差并加上均值，得到原始数据空间上的PCA之后数据。下图展示了这五步进行主成分分析的过程。

2.5 主成分分析的python应用

在python的经典机器学习库中，主成分分析的模型函数对应sklearn.decomposition.PCA函数，具体为sklearn.decomposition.PCA(n_components=None, *, copy=True, whiten=False, svd_solver='auto', tol=0.0, iterated_power='auto', random_state=None)。其中n_components是要保留的组件数，即主成分的个数，默认保留所有组件；copy值如果为False，则传递给fit的数据将被覆盖，并且运行fit(X).transform(X)将不会产生预期的结果，请改用fit_transform(X)；whiten值如果为True，则将component_向量乘以n_samples的平方根，然后除以奇异值，以确保具有单位分量方差的不相关输出；svd_solver是设置奇异值分解（特征分解）的方法，如果是默认的auto则通过基于X.shape和n_components的默认策略选择求解器，即如果输入数据大于500x500并且要提取的组件数小于数据最小维度的80％，则选择“随机化”方法，否则将计算出精确的完整SVD（SVD是Singular Value Decomposition，奇异值分解）再选择截断。下面是一个PCA的python例子。

import numpy as np
from sklearn.decomposition import PCA
X = np.array([[-1, -1], [-2, -1], [-3, -2], [1, 1], [2, 1], [3, 2]])  #原始数据
pca = PCA(n_components=2)  #选择主成分个数数为2的PCA模型
pca.fit(X)  #向PCA模型输入数据
print(pca.explained_variance_ratio_) #打印出由每个主成分占有的方差百分比。
#结果是[0.9924... 0.0075...]，由高到低排列
print(pca.singular_values_) #打印出每个主成分的特征值
#结果是[6.30061... 0.54980...]

#重新设置PCA模型的参数
pca = PCA(n_components=2, svd_solver='full')  #将PCA模型中的SVD参数设置为full，即通过scipy.linalg.svd运行精确的完整SVD调用标准LAPACK求解器
pca.fit(X)
print(pca.explained_variance_ratio_)
#百分比结果是[0.9924... 0.00755...]
print(pca.singular_values_)
#特征值结果是[6.30061... 0.54980...]

python的sklearn库也提供了基于PCA扩展之后的方法，比如KernelPCA，它是PCA的扩展，它通过使用内核来实现非线性降维，可以用于降噪，压缩和结构化预测（内核相关性估计）。KernelPCA的函数是sklearn.decomposition.KernelPCA()，具体可见KernelPCA官方手册，这里不做过多解释。

3. 聚类（Clustering）

3.1 聚类分析的简介

聚类分析（cluster analysis 或 clustering）是一种非监督性学习，原始数据是不带标签，即不含有人为添加的信息。聚类分析将原始数据点根据距离或者相似度将其分为几个特定的组或者聚类（cluster）中，结果使得同一组中的对象彼此之间（在某种意义上）比其他组中的对象更相似，与其他组有明显的不同。聚类分析本身不是一种特定的算法，而是要解决的一般任务，而聚类公式可以化为一个多目标优化问题。其过程可以通过各种算法来实现，常见算法包括kmeans算法，DBSCAN聚类算法，高斯混合模型，和层次聚类算法。

聚类分析主要分为聚集型聚类（Agglomerative）和分散性聚类（Divisive）。聚集型聚类是将成对的数据点（items）或者聚类相继连接以产生更大的聚类。分散性聚类是将所有数据点最初放置在一个群集中，然后分成单独的聚类。聚集型聚类的代表算法是分层聚类（Hierarchical clustering），将在3.3节中详细讨论。分散性聚类的代表算法是K均值（Kmeans）算法，将在3.4节中详细讨论。

3.2 聚类中的距离和相似度

除了各种具体算法，聚类分析最重要的一步是选择要测量的距离，或者是数据点之间的相似度测量标准，以此距离或标准决定之后划归具体聚类。最简单的测量距离方法是使用曼哈顿距离（Manhattan distance），即每个变量的绝对差值之和；曼哈顿距离也叫L1距离，1-norm距离，其本质和平均绝对误差损失函数（L1 Loss）相同，具体L1损失函数参见上一篇文章损失函数一节。用的最多的距离是欧几里德距离（Euclidean Distance，简称欧式距离），也叫2-norm距离，使用的是欧几里得空间中两点间“普通”（即直线）距离，它的数学计算公式为

计算欧式距离要找到一个空间，来计算每个空间中数据点到原点的距离，然后对所有距离进行换算。曼哈顿距离或欧式距离越小的点相距约近，更有可能被分到同一个聚类中。除了使用距离，还可以使用相似度判断每个点的归属，常用的是余弦相似度（Euclidean Distance）。余弦相似度通过测量两个向量的夹角的余弦值来度量它们之间的相似性。其中0度角的余弦值是1，而其他任何角度的余弦值都不大于1，并且不小于-1。两个向量之间的角度的余弦值确定两个向量是否大致指向相同的方向，所以余弦相似度越大，数据点越相关，更有可能属于同一个聚类。余弦相似度的数学计算公式为

通过计算每个数据点之间的距离或相似度，可以得到所有点的相似度（距离）矩阵，其中矩阵中的第（i，j）项是数据点i和j之间的距离（相似度），而且这个矩阵是对角线对称矩阵，相似度矩阵的对角线数均为1，距离矩阵的对角线数均为0，所以只需要求解出这个相似度矩阵的下三角矩阵。相似度矩阵的一个实例是文档的术语相似度矩阵（Term Similarities Matrix），多用于自然语言处理中的文档处理和数据检索。下面以一个文档分析中的术语相似度为例，展示相似度（距离）进行矩阵计算和图像化的过程。下面矩阵是各个术语（term）在每个文档（document）中出现的频率

术语相似度的计算公式为

根据这个公式计算得到上述表格的术语相似度矩阵如下，其中仅需计算矩阵的下三角部分，相似度越大表明这两个术语之间越相似，包含更多的相同信息，根据这个矩阵可以分析每个术语的重要程度，便于进一步计算tf-idf值。

上述矩阵虽然可以具体显示出术语之间的相似程度，但是不便于直接观察哪几个术语是相似的。这是可以设置一个相似度（或距离）阈值（threshold），当两个数据点（此处为术语term）的相似度大于这个阈值则判断这两个数据点是相关或者近似的，将其记为1，否则记为0。那么上面的术语相似度矩阵在阈值设为10的情况下可以化为下面矩阵，其中的1值可以直接看出哪两个术语之间的相关的。

为了视觉上更加直观，相似度矩阵可以进行图像展示（Graph Representation），相似度矩阵通过一个无向图（undirected graph）进行视觉化，其中每一个数据（这里是术语）化为一个图像上的点（node），如果经过阈值处理后的相似度矩阵中这两个数据点的值为1，则判断这两个点是相似的，并在图中用一个没有箭头方向的直线将这两个点相连接。如果相似度矩阵没有阈值，可以在图像中加入权重，代表这两个点的相似程度。上面实例的术语相似度矩阵图像化后变为

3.3 分层聚类

分层聚类（Hierarchical clustering）也叫Hierarchical Agglomerative方法，是通过依次合并或拆分嵌套的聚类来构建最终的聚类。这种群集的层次结构可以表示为树（或树状图），树的根是收集所有样本的唯一簇，叶子是只有一个样本的簇。分层聚类的树状图表示如下图

分层聚类算法的计算过程是：首先计算距离矩阵，默认使用欧式距离，也可以使用相似度矩阵；然后将每个数据点放在自己的群集中；接着查找并合并最相似的集群对，以树状图显示合并；之后更新相似度（距离）矩阵。重复上述步骤直到所有数据点都在一个群集中。

分层聚类在python中的对应函数为sklearn.cluster.AgglomerativeClusterin。具体为sklearn.cluster.AgglomerativeClustering(n_clusters=2, *, affinity='euclidean', memory=None, connectivity=None, compute_full_tree='auto', linkage='ward', distance_threshold=None, compute_distances=False)，其中n_clusters是要查找的集群（聚类）数，默认为2；affinity用于计算距离或者相似度的度量，可以选择 “euclidean”, “l1”, “l2”, “manhattan”, “cosine”；memory用于缓存树计算的输出；compute_full_tree设置尽早在n_clusters处停止树的构建，如果簇的数量与样本的数量相比不小的话，这可以减少计算时间很；linkage代表使用哪个链接标准，用于确定要在观察组之间使用的距离，最后将合并最小化此标准的成对聚类。下面是一个简单实例。

from sklearn.cluster import AgglomerativeClustering
import numpy as np
X = np.array([[1, 2], [1, 4], [1, 0],[4, 2], [4, 4], [4, 0]]) #原始数据点
clustering = AgglomerativeClustering().fit(X) #输入数据
clustering.labels_   #打印每个点的聚类标签
#聚类的结果显示每个点的所属标签，array([1, 1, 1, 0, 0, 0])

3.4 K均值算法

K均值（KMeans）算法以空间中k个点为中心进行聚类，通过将最靠近中心的数据点聚合，把所有点分成等方差的k个集群，每个集群内平方和都最小化。KMeans算法简单，计算快速，适用于大样本数据。每个集群的平均点称为中心（centres）或质点（centroids），使用μ表示，但是质点通常不是来自原始数据集X，虽然它们在同一空间。K-means算法旨在选择最小化惯性（inertia）的中心，和群集内平方和（sum-of-squares）的标准。所以K-means算法的核心步骤是重复下面两步直到收敛：对于每个集群找到最合适的中心，对于每个中心找到符合标准的属于它的集群点。

KMeans算法（也称Lloyd’s algorithm）具体步骤是：选择聚类的个数k；任意产生k个聚类，然后确定聚类中心，或者直接生成k个中心（可以随机生成）；对每个点确定其所属的聚类中心点；根据新的聚类情况再计算其聚类新中心。重复以上步骤直到满足收敛要求，通常是确定的中心点不再改变。初始化中随机中心点的生成对结果有着重大影响，所以kmeans算法不稳定，聚类结果可能比较差，不能够满足需要。

KMeans算法等效于带有最小，全等，对角协方差矩阵的期望最大化算法（EM算法），EM算法将在下面4.3节中讲解。按照EM算法步骤求解kmeans，先确定需要最小化的损失函数，kmeans中要使得群集内平方和最小，所以损失函数是所有集群中点到中心的距离平方之和，即

EM算法中对应的初始化步骤为初始化K个中心点μ₁ ,…, μ_k 。E步为更新集群，根据距离数据点最近的中心点判断所属集群，使用r_nk 来表示每个点属于哪个集群，如果属于则为1，否则为0，数学计算式为

M步为更新中心点，根据新的集群计算新的平均值点作为新的中心点，对于k个中心点都通过下面式子更新

重复E步和M步直到中心点收敛，或者中心点移动的距离小于阈值（相对公差）时迭代停止。

训练有足够的时间，K均值算法将始终收敛，但是这可能是收敛到局部最小值。是否收敛到局部最小值在很大程度上取决于中心点的初始化，多次随机初始化中心点会改善结果。另一个办法是使用k-means ++初始化方案，将初始化的中心点彼此相离尽可能远，效果好于随机初始化，也可以单独调用K-means ++来选择其他聚类算法的种子作为初始化的中心点。K均值算法的结果也受超参数k的影响，可以通过比较训练集数据和验证集的损失函数结果，来选择合适的k值，一般是在训练集数据的损失函数下降曲线的拐点选择k。

K均值算法在python中的对应函数为sklearn.cluster.KMeans()。具体为sklearn.cluster.KMeans(n_clusters=8, *, init='k-means++', n_init=10, max_iter=300, tol=0.0001, precompute_distances='deprecated', verbose=0, random_state=None, copy_x=True, n_jobs='deprecated', algorithm='auto')，其中n_clusters为集群的数量以及要生成的质心的数量，默认值为8；init是调用初始化的方法，默认为k-means ++，也可以选择随机生成；n_init代表k均值算法将在不同质心种子下运行的时间；max_iter是单次运行的k均值算法的最大迭代次数；tol是关于Frobenius范数的相对公差，代表中心收敛判断的阈值；random_state用于质心初始化的随机数生成。具体实例为

from sklearn.cluster import KMeans
import numpy as np
X = np.array([[1, 2], [1, 4], [1, 0],[10, 2], [10, 4], [10, 0]])
kmeans = KMeans(n_clusters=2, random_state=0).fit(X)
kmeans.labels_
#结果是array([1, 1, 1, 0, 0, 0], dtype=int32)
kmeans.predict([[0, 0], [12, 3]])
#结果是array([1, 0], dtype=int32)
kmeans.cluster_centers_
#结果是array([[10.,  2.],[ 1.,  2.]])

4.高斯混合模型（Gaussian Mixture Models）

4.1 高斯混合模型（GMM）的由来

高斯分布是最常见的连续随机变量的分布，自然和社会中很多的概率分布都是高斯分布，或者是几个高斯分布组成。使用多个高斯分布可以近似表达出几乎任意分布状态，即可以平滑任意形状的概率分布。高斯混合模型（Gaussian Mixture Models，简称GMM）就是一个随机概率分布分解成几个高斯分布的联合表示过程，找到这些高斯分布线性组合的方式。可以使用高斯混合模型进行密度估计。而理论上高斯混合模型可以估计表述出任意模型。以下图所示，下图的分布可以看成由红，蓝和绿三个不同分布相加而成，每一个分布都可以近似看作是一个简单的分布。

总体分布的概率密度等于各个单一分布密度乘以各自的加权值后相加，数学表达式如下。其中p(x)是要描述的分布的概率密度，它由K个单一分布组成，p_k是每一个单一分布（basic distributions）的概率（可以是高斯分布Gaussians，伯努利分布Bernoullis或伽马分布Gammas），π_k是混合的权重，权重和为1。

下图展示了一个由三个高斯分布相加得到的一个混合分布，包含蓝色，黄色和绿色三个高斯分布，而黑色曲线是混合分布的概率曲线。

这个高斯混合分布的概率密度函数可以使用三个单一高斯分布的加权相加得到，数学表达式为（其中不同颜色部分对应图中相同颜色的分布）

而高斯混合模型就是求解上述式子的加权值和各个高斯分布的参数（均值和标准差），即求解下面式子中的π_k，μ_k和Σ_k，这些参数统称为θ，也是高斯混合模型的参数集合。

高斯混合模型对比单一高斯分布，它有更大的弹性，适用更复杂的分布变化，可以处理更复杂的模型分布。求解高斯混合模型的目的是最大化表述整体分布，也就是高斯混合模型可以最大化拟合分布曲线，这可以参考线性回归中为了使得拟合最大化所采用的最大似然方法。事实上，高斯混合模型在数学上就是使用最大似然方法进行求解。

4.2 最大似然估计求解GMM

对于数据集X={x₁, … , x_N}，根据高斯混合模型，总体的概率密度函数等于各个高斯概率函数之和，即（其中p(X)是总的概率函数，p(x)是单一高斯概率）

下面使用最大似然估计求解模型的系数θ。首先似然概率是上面式子所示的概率乘积，对这个似然概率求对数之后得到对数似然概率（log-likelihood），也就是求对数之后的损失函数

为了使得高斯混合模型最佳表现出整体概率分布，损失函数要最小，那么求对数之后的损失函数则是要最大，所以要找到π_k，μ_k和Σ_k组合使得上述对数似然概率最大。将对数似然概率分别对π_k，μ_k和Σ_k求偏导数并使其为0，即可以达到对数似然概率最大。

根据链式法则将三个偏导数式子混合成一个参数θ={π_k，μ_k和Σ_k: k=1, …, K}，求导式子变为

其中黄色部分为

为了式子表达简介，定义一个术语叫quantity，也叫responsibilities，使用r_nk表示，它可以看作第n个数据点的第k个高斯混合分量的所占比例。具体为

由r_nk定义可以得到每个r_nk到大于或等于0，而且同一个数据点（n一样）的总共K个r_nk之和为1。而根据刚才导数为0的式子可以把模型系数π_k，μ_k和Σ_k的解使用r_nk表示，为（具体数学过程省略）

4.3 EM算法

使用最大似然估计可以得到高斯混合模型参数的数学求解，但是在实际应用中使用EM算法在迭代过程中求解参数。EM（Expectation-Maximum）算法也称期望最大化算法，常被用来学习高斯混合模型的参数，隐式马尔科夫算法（HMM）或LDA主题模型的变分推断等。EM算法是一种迭代优化策略，每一次迭代分为两步，首先是期望步（E步，Expectation），然后是极大步（M步，Maximum）。EM算法和极大似然估计的使用前提是一样的，需要假设数据总体的分布，如果不知道数据分布则无法使用EM算法。所以EM算法也被看作一个在已知部分相关变量的情况下，估计未知变量的迭代技术。而EM算法在迭代之前就需要初始化分布参数。

EM算法的步骤使用流程图表示为

具体步骤为：首先初始化一组参数的值。 假设观测数据来自特定模型，则会向系统提供一组不完整的观测数据。然后进行E步（E-step），使用观察到的数据来估计或猜测缺失或不完整数据的值，进行更新变量。接下来是M步（M-step），使用在前面的E步中生成的完整数据，以更新参数的值，同时更新假设。最后检查参数值是否收敛，如果是则停止迭代，否则重复步骤2和步骤3，即E步和M步，直到收敛为止。

在求解高斯混合模型的过程中，EM算法计算步骤是：首先初始化高斯混合模型参数π_k，μ_k和Σ_k，其中一种初始化方法为π_k=1/K，μ_k取k均值算法后的中心centroids，Σ_k取方差值。然后进行E步，使用3.2节中倒数第二个式子即r_nk的式子计算权重（responsibilities）r_nk的值，也就是属于混合分量k的数据点n的后验概率。接下来进行M步，使用估计得到的r_nk更新模型参数π_k，μ_k和Σ_k，即使用3.2节中最后的π_k，μ_k和Σ_k的三个解计算式子更新参数。最后重复E步和M步直到参数的值收敛。

高斯混合模型求解最后可以收敛于局部最小值，但是因为它使用软标签（soft labels），即同一个数据点可以属于多个高斯分布，所以高斯混合模型的局部最小值也处于较优状态，软标签允许数据点在不同簇（clusters）之间缓慢移动。

4.4 高斯混合模型的优缺点

高斯混合模型的优点是：GMM算法的速度快，是学习混合模型的最快算法；而且GMM算法具有Agnostic性，由它仅使可能性最大化，因此不会使均值趋于零，也不会使聚类大小具有可能适用或可能不适用的特定结构。

高斯混合模型的缺点是：GMM算法具有奇点（Singularities），当每个混合点的点数不足时，估计协方差矩阵将变得困难，并且除非有人为地对协方差进行人为调整，否则该算法会发散并找到具有无限可能性的解；而且GMM算法的组件数（Number of components）可能过大，这个算法将始终使用它可以访问的所有组件（各个单一高斯分布），需要保留的数据或信息理论标准来确定在没有外部提示的情况下要使用多少个组件，导致最终组件数过大。

4.5 高斯混合模型的python应用

高斯混合模型在python的sklearn库中对应sklearn.mixture.GaussianMixture()函数。具体为sklearn.mixture.GaussianMixture(n_components=1, *, covariance_type='full', tol=0.001, reg_covar=1e-06, max_iter=100, n_init=1, init_params='kmeans', weights_init=None, means_init=None, precisions_init=None, random_state=None, warm_start=False, verbose=0, verbose_interval=10)，其中n_components是高斯混合成分的数量，默认为1；covariance_type为描述要使用的协方差参数类型的字符串，可以是‘full’， ‘tied’， ‘diag’或‘spherical’；tol是收敛阈值；reg_covar代表非负正则化项添加到协方差矩阵的对角线上，确保协方差矩阵均为正；max_iter是EM算法最大迭代次数；n_init是进行初始化的次数，从中选择最佳结果作为最终初始化；init_params是选择用于初始化权重，均值和精度的方法，如果是‘kmeans’则使用kmeans方法初始化，如果是’random’则是随机初始化；weights_init是用户提供的初始权重，如果为None则使用init_params方法生成的；means_init和precisions_init分别是用户提供的初始均值，和初始精度（协方差矩阵的倒数）；random_state是控制分配给用于初始化参数的方法的随机种子；warm_start默认为False，如果为True则将最后一次拟合的结果用作对fit()的下一次调用的初始化。下面是一个实例。

import numpy as np
from sklearn.mixture import GaussianMixture
X = np.array([[1, 2], [1, 4], [1, 0], [10, 2], [10, 4], [10, 0]]) #原始数据
gm = GaussianMixture(n_components=2, random_state=0).fit(X) #选择成分数为2的高斯混合模型，使用0号随机种子初始化，将数据X输入模型
gm.means_ #输出每种高斯混合成分的平均值
#均值结果是  array([[10.,  2.],
#                 [ 1.,  2.]])
#使用gm.covariances_ 可以输出每种高斯混合成分的方差
#使用gm.precisions_ 可以输出高斯混合模型中每种成分的精密度矩阵
gm.predict([[0, 0], [12, 3]])  #输入新的数据进行预测，生成数据点所属高斯分布的预测标签
#预测结果是 array([1, 0])