用Python实现降维和聚类

主成分分析（Principal components analysis）-最大方差解释

     在这一篇之前的内容是《Factor Analysis》，由于非常理论，打算学完整个课程后再写。在写这篇之前，我阅读了PCA、SVD和LDA。这几个模型相近，却都有自己的特点。本篇打算先介绍PCA，至于他们之间的关系，只能是边学边体会了。PCA以前也叫做Principal factor analysis。

1. 问题

     真实的训练数据总是存在各种各样的问题：

1、 比如拿到一个汽车的样本，里面既有以“千米/每小时”度量的最大速度特征，也有“英里/小时”的最大速度特征，显然这两个特征有一个多余。

2、 拿到一个数学系的本科生期末考试成绩单，里面有三列，一列是对数学的兴趣程度，一列是复习时间，还有一列是考试成绩。我们知道要学好数学，需要有浓厚的兴趣，所以第二项与第一项强相关，第三项和第二项也是强相关。那是不是可以合并第一项和第二项呢？

3、 拿到一个样本，特征非常多，而样例特别少，这样用回归去直接拟合非常困难，容易过度拟合。比如北京的房价：假设房子的特征是（大小、位置、朝向、是否学区房、建造年代、是否二手、层数、所在层数），搞了这么多特征，结果只有不到十个房子的样例。要拟合房子特征->房价的这么多特征，就会造成过度拟合。

4、 这个与第二个有点类似，假设在IR中我们建立的文档-词项矩阵中，有两个词项为“learn”和“study”，在传统的向量空间模型中，认为两者独立。然而从语义的角度来讲，两者是相似的，而且两者出现频率也类似，是不是可以合成为一个特征呢？

5、 在信号传输过程中，由于信道不是理想的，信道另一端收到的信号会有噪音扰动，那么怎么滤去这些噪音呢？

     回顾我们之前介绍的《模型选择和规则化》，里面谈到的特征选择的问题。但在那篇中要剔除的特征主要是和类标签无关的特征。比如“学生的名字”就和他的“成绩”无关，使用的是互信息的方法。

     而这里的特征很多是和类标签有关的，但里面存在噪声或者冗余。在这种情况下，需要一种特征降维的方法来减少特征数，减少噪音和冗余，减少过度拟合的可能性。

     下面探讨一种称作主成分分析（PCA）的方法来解决部分上述问题。PCA的思想是将n维特征映射到k维上（k<n），这k维是全新的正交特征。这k维特征称为主元，是重新构造出来的k维特征，而不是简单地从n维特征中去除其余n-k维特征。

2. PCA计算过程

     首先介绍PCA的计算过程：

     假设我们得到的2维数据如下：

     

     行代表了样例，列代表特征，这里有10个样例，每个样例两个特征。可以这样认为，有10篇文档，x是10篇文档中“learn”出现的TF-IDF，y是10篇文档中“study”出现的TF-IDF。也可以认为有10辆汽车，x是千米/小时的速度，y是英里/小时的速度，等等。

     第一步分别求x和y的平均值，然后对于所有的样例，都减去对应的均值。这里x的均值是1.81，y的均值是1.91，那么一个样例减去均值后即为（0.69,0.49），得到

     

     第二步，求特征协方差矩阵，如果数据是3维，那么协方差矩阵是

     

     这里只有x和y，求解得

     

     对角线上分别是x和y的方差，非对角线上是协方差。协方差大于0表示x和y若有一个增，另一个也增；小于0表示一个增，一个减；协方差为0时，两者独立。协方差绝对值越大，两者对彼此的影响越大，反之越小。

     第三步，求协方差的特征值和特征向量，得到

     

     上面是两个特征值，下面是对应的特征向量，特征值0.0490833989对应特征向量为，这里的特征向量都归一化为单位向量。

    第四步，将特征值按照从大到小的顺序排序，选择其中最大的k个，然后将其对应的k个特征向量分别作为列向量组成特征向量矩阵。

     这里特征值只有两个，我们选择其中最大的那个，这里是1.28402771，对应的特征向量是。

     第五步，将样本点投影到选取的特征向量上。假设样例数为m，特征数为n，减去均值后的样本矩阵为DataAdjust(m*n)，协方差矩阵是n*n，选取的k个特征向量组成的矩阵为EigenVectors(n*k)。那么投影后的数据FinalData为

     

     这里是

     FinalData(10*1) = DataAdjust(10*2矩阵)×特征向量

     得到结果是

     

     这样，就将原始样例的n维特征变成了k维，这k维就是原始特征在k维上的投影。

     上面的数据可以认为是learn和study特征融合为一个新的特征叫做LS特征，该特征基本上代表了这两个特征。

     上述过程有个图描述：

     

     正号表示预处理后的样本点，斜着的两条线就分别是正交的特征向量（由于协方差矩阵是对称的，因此其特征向量正交），最后一步的矩阵乘法就是将原始样本点分别往特征向量对应的轴上做投影。

     如果取的k=2，那么结果是

     

     这就是经过PCA处理后的样本数据，水平轴（上面举例为LS特征）基本上可以代表全部样本点。整个过程看起来就像将坐标系做了旋转，当然二维可以图形化表示，高维就不行了。上面的如果k=1，那么只会留下这里的水平轴，轴上是所有点在该轴的投影。

     这样PCA的过程基本结束。在第一步减均值之后，其实应该还有一步对特征做方差归一化。比如一个特征是汽车速度（0到100），一个是汽车的座位数（2到6），显然第二个的方差比第一个小。因此，如果样本特征中存在这种情况，那么在第一步之后，求每个特征的标准差，然后对每个样例在该特征下的数据除以。

     归纳一下，使用我们之前熟悉的表示方法，在求协方差之前的步骤是：

     

     其中是样例，共m个，每个样例n个特征，也就是说是n维向量。是第i个样例的第j个特征。是样例均值。是第j个特征的标准差。

     整个PCA过程貌似及其简单，就是求协方差的特征值和特征向量，然后做数据转换。但是有没有觉得很神奇，为什么求协方差的特征向量就是最理想的k维向量？其背后隐藏的意义是什么？整个PCA的意义是什么？

3. PCA理论基础

     要解释为什么协方差矩阵的特征向量就是k维理想特征，我看到的有三个理论：分别是最大方差理论、最小错误理论和坐标轴相关度理论。这里简单探讨前两种，最后一种在讨论PCA意义时简单概述。

3.1 最大方差理论

     在信号处理中认为信号具有较大的方差，噪声有较小的方差，信噪比就是信号与噪声的方差比，越大越好。如前面的图，样本在横轴上的投影方差较大，在纵轴上的投影方差较小，那么认为纵轴上的投影是由噪声引起的。

因此我们认为，最好的k维特征是将n维样本点转换为k维后，每一维上的样本方差都很大。

     比如下图有5个样本点：（已经做过预处理，均值为0，特征方差归一）

     

     下面将样本投影到某一维上，这里用一条过原点的直线表示（前处理的过程实质是将原点移到样本点的中心点）。

     

     假设我们选择两条不同的直线做投影，那么左右两条中哪个好呢？根据我们之前的方差最大化理论，左边的好，因为投影后的样本点之间方差最大。

     这里先解释一下投影的概念：

     

     红色点表示样例，蓝色点表示在u上的投影，u是直线的斜率也是直线的方向向量，而且是单位向量。蓝色点是在u上的投影点，离原点的距离是（即或者）由于这些样本点（样例）的每一维特征均值都为0，因此投影到u上的样本点（只有一个到原点的距离值）的均值仍然是0。

     回到上面左右图中的左图，我们要求的是最佳的u，使得投影后的样本点方差最大。

     由于投影后均值为0，因此方差为：

     

     中间那部分很熟悉啊，不就是样本特征的协方差矩阵么（的均值为0，一般协方差矩阵都除以m-1，这里用m）。

     用来表示，表示，那么上式写作

      

     由于u是单位向量，即，上式两边都左乘u得，

     即

     We got it！就是的特征值，u是特征向量。最佳的投影直线是特征值最大时对应的特征向量，其次是第二大对应的特征向量，依次类推。

     因此，我们只需要对协方差矩阵进行特征值分解，得到的前k大特征值对应的特征向量就是最佳的k维新特征，而且这k维新特征是正交的。得到前k个u以后，样例通过以下变换可以得到新的样本。

     

     其中的第j维就是在上的投影。

     通过选取最大的k个u，使得方差较小的特征（如噪声）被丢弃。

3.2 最小平方误差理论

     

     假设有这样的二维样本点（红色点），回顾我们前面探讨的是求一条直线，使得样本点投影到直线上的点的方差最大。本质是求直线，那么度量直线求的好不好，不仅仅只有方差最大化的方法。再回想我们最开始学习的线性回归等，目的也是求一个线性函数使得直线能够最佳拟合样本点，那么我们能不能认为最佳的直线就是回归后的直线呢？回归时我们的最小二乘法度量的是样本点到直线的坐标轴距离。比如这个问题中，特征是x，类标签是y。回归时最小二乘法度量的是距离d。如果使用回归方法来度量最佳直线，那么就是直接在原始样本上做回归了，跟特征选择就没什么关系了。

     因此，我们打算选用另外一种评价直线好坏的方法，使用点到直线的距离d’来度量。

     现在有n个样本点，每个样本点为m维（这节内容中使用的符号与上面的不太一致，需要重新理解符号的意义）。将样本点在直线上的投影记为，那么我们就是要最小化

     

     这个公式称作最小平方误差（Least Squared Error）。

     而确定一条直线，一般只需要确定一个点，并且确定方向即可。

     第一步确定点：

     假设要在空间中找一点来代表这n个样本点，“代表”这个词不是量化的，因此要量化的话，我们就是要找一个m维的点，使得

     

     最小。其中是平方错误评价函数（squared-error criterion function），假设m为n个样本点的均值：

     

     那么平方错误可以写作：

     

     后项与无关，看做常量，而，因此最小化时，

      

     是样本点均值。

     第二步确定方向：

     我们从拉出要求的直线（这条直线要过点m），假设直线的方向是单位向量e。那么直线上任意一点，比如就可以用点m和e来表示

      

     其中是到点m的距离。

     我们重新定义最小平方误差：

     

     这里的k只是相当于i。就是最小平方误差函数，其中的未知参数是和e。

     实际上是求的最小值。首先将上式展开：

     

     我们首先固定e，将其看做是常量，，然后对进行求导，得

     

     这个结果意思是说，如果知道了e，那么将与e做内积，就可以知道了在e上的投影离m的长度距离，不过这个结果不用求都知道。

     然后是固定，对e求偏导数，我们先将公式（8）代入，得 

     

     其中 与协方差矩阵类似，只是缺少个分母n-1，我们称之为散列矩阵（scatter matrix）。

     然后可以对e求偏导数，但是e需要首先满足，引入拉格朗日乘子，来使最大（最小），令

     

     求偏导

     

     这里存在对向量求导数的技巧，方法这里不多做介绍。可以去看一些关于矩阵微积分的资料，这里求导时可以将看作是，将看做是。

     导数等于0时，得

     

     两边除以n-1就变成了，对协方差矩阵求特征值向量了。

     从不同的思路出发，最后得到同一个结果，对协方差矩阵求特征向量，求得后特征向量上就成为了新的坐标，如下图：

     

     这时候点都聚集在新的坐标轴周围，因为我们使用的最小平方误差的意义就在此。

4. PCA理论意义

     PCA将n个特征降维到k个，可以用来进行数据压缩，如果100维的向量最后可以用10维来表示，那么压缩率为90%。同样图像处理领域的KL变换使用PCA做图像压缩。但PCA要保证降维后，还要保证数据的特性损失最小。再看回顾一下PCA的效果。经过PCA处理后，二维数据投影到一维上可以有以下几种情况：

     

     我们认为左图好，一方面是投影后方差最大，一方面是点到直线的距离平方和最小，而且直线过样本点的中心点。为什么右边的投影效果比较差？直觉是因为坐标轴之间相关，以至于去掉一个坐标轴，就会使得坐标点无法被单独一个坐标轴确定。

     PCA得到的k个坐标轴实际上是k个特征向量，由于协方差矩阵对称，因此k个特征向量正交。看下面的计算过程。

     假设我们还是用来表示样例，m个样例，n个特征。特征向量为e，表示第i个特征向量的第1维。那么原始样本特征方程可以用下面式子来表示：

     前面两个矩阵乘积就是协方差矩阵（除以m后），原始的样本矩阵A是第二个矩阵m*n。

     

     上式可以简写为

     我们最后得到的投影结果是，E是k个特征向量组成的矩阵，展开如下：

     

     得到的新的样例矩阵就是m个样例到k个特征向量的投影，也是这k个特征向量的线性组合。e之间是正交的。从矩阵乘法中可以看出，PCA所做的变换是将原始样本点（n维），投影到k个正交的坐标系中去，丢弃其他维度的信息。举个例子，假设宇宙是n维的（霍金说是11维的），我们得到银河系中每个星星的坐标（相对于银河系中心的n维向量），然而我们想用二维坐标去逼近这些样本点，假设算出来的协方差矩阵的特征向量分别是图中的水平和竖直方向，那么我们建议以银河系中心为原点的x和y坐标轴，所有的星星都投影到x和y上，得到下面的图片。然而我们丢弃了每个星星离我们的远近距离等信息。

     

5. 总结与讨论

     这一部分来自http://www.cad.zju.edu.cn/home/chenlu/pca.htm

     PCA技术的一大好处是对数据进行降维的处理。我们可以对新求出的“主元”向量的重要性进行排序，根据需要取前面最重要的部分，将后面的维数省去，可以达到降维从而简化模型或是对数据进行压缩的效果。同时最大程度的保持了原有数据的信息。

     PCA技术的一个很大的优点是，它是完全无参数限制的。在PCA的计算过程中完全不需要人为的设定参数或是根据任何经验模型对计算进行干预，最后的结果只与数据相关，与用户是独立的。 
但是，这一点同时也可以看作是缺点。如果用户对观测对象有一定的先验知识，掌握了数据的一些特征，却无法通过参数化等方法对处理过程进行干预，可能会得不到预期的效果，效率也不高。

     

     图表 4：黑色点表示采样数据，排列成转盘的形状。 
     容易想象，该数据的主元是或是旋转角。

     如图表 4中的例子，PCA找出的主元将是。但是这显然不是最优和最简化的主元。之间存在着非线性的关系。根据先验的知识可知旋转角是最优的主元（类比极坐标）。则在这种情况下，PCA就会失效。但是，如果加入先验的知识，对数据进行某种划归，就可以将数据转化为以为线性的空间中。这类根据先验知识对数据预先进行非线性转换的方法就成为kernel-PCA，它扩展了PCA能够处理的问题的范围，又可以结合一些先验约束，是比较流行的方法。

     有时数据的分布并不是满足高斯分布。如图表 5所示，在非高斯分布的情况下，PCA方法得出的主元可能并不是最优的。在寻找主元时不能将方差作为衡量重要性的标准。要根据数据的分布情况选择合适的描述完全分布的变量，然后根据概率分布式

     

     来计算两个向量上数据分布的相关性。等价的，保持主元间的正交假设，寻找的主元同样要使。这一类方法被称为独立主元分解(ICA)。

     

     图表 5：数据的分布并不满足高斯分布，呈明显的十字星状。 
     这种情况下，方差最大的方向并不是最优主元方向。

     另外PCA还可以用于预测矩阵中缺失的元素。

6. 其他参考文献

     A tutorial on Principal Components Analysis LI Smith – 2002

     A Tutorial on Principal Component Analysis J Shlens

     http://www.cmlab.csie.ntu.edu.tw/~cyy/learning/tutorials/PCAMissingData.pdf

     http://www.cad.zju.edu.cn/home/chenlu/pca.htm

K-Means算法

非监督式学习对一组无标签的数据试图发现其内在的结构，主要用途包括：

• 市场划分（Market Segmentation）
• 社交网络分析（Social Network Analysis）
• 管理计算机集群（Organize Computer Clusters）
• 天文学数据分析（Astronomical Data Analysis）

K-Means算法属于非监督式学习的一种，算法的输入是：训练数据集 {x(1),x(2),…,x(m)} {x(1),x(2),…,x(m)}(其中 x(i)∈Rn x(i)∈Rn)和聚类数量 K K（将数据划分为 K K类）；算法输出是 K K个聚类中心 μ1,μ2,…,μK μ1,μ2,…,μK和每个数据点 x(i) x(i)所在的分类。

K-Means算法步骤

1. 随机初始化 K K个聚类中心(cluster centroid)  μ1,μ2,…,μK μ1,μ2,…,μK
2. Cluster Assignment: 对于每个数据点 x(i) x(i)，寻找离它最近的聚类中心，将其归入该类；即 c(i)=mink||x(i)−μk||2 c(i)=mink||x(i)−μk||2，其中 c(i) c(i)表示 x(i) x(i)所在的类
3. Move Centroid: 更新聚类中心 uk uk的值为所有属于类 k k的数据点的平均值
4. 重复2、3步直到收敛或者达到最大迭代次数

图1 K-Means算法示例

K-Means算法的优化目标

用 μc(i) μc(i)表示第 i i个数据点 x(i) x(i)所在类的中心，则K-Means优化的代价函数为

J(c(1),…,c(m),μ1,…,μK)=1m∑i=1m||x(i)−μc(i)||2 J(c(1),…,c(m),μ1,…,μK)=1m∑i=1m||x(i)−μc(i)||2

希望找到最优参数使得该函数最小化，即

minc(1),…,c(m)μ1,…,μKJ(c(1),…,c(m),μ1,…,μK) minc(1),…,c(m)μ1,…,μKJ(c(1),…,c(m),μ1,…,μK)

需要注意的问题

• 随机初始化：常用的初始化方法是，从训练数据点中随机选择 K K( K<m K<m)个数据点，作为初始的聚类中心 μ1,μ2,…,μK μ1,μ2,…,μK
• 局部最优：算法聚类的性能与初始聚类中心的选择有关，为避免陷入局部最优（如图2所示），应该运行多次（50次）取使得 J J最小的结果
• K K值选择：Elbow方法，绘制 J J随 K K的变化曲线，选择下降速度突然变慢的转折点作为K值；对于转折不明显的曲线，可根据K-Means算法后续的目标进行选择。

   
图2 K-Means算法的全局最优解和局部最优解    
                  

图3  用Elbow方法选择K值的情况（左）和Elbow法不适用的情况（右）

 PCA降维算法

动机

数据压缩：将高维数据(n维)压缩为低维数据(k维)

数据可视化：将数据压缩到2维/3维方便可视化

PCA问题形式化

如果需要将二维数据点，压缩为一维数据点，我们需要找到一个方向，使得数据点到这个方向上投射时的误差最小（即点到该直线的距离最小）；更一般地，如果需要将 n n维的数据点压缩到 k k维，我们需要找到 k k个新的方向 u(1),u(2),…,u(k) u(1),u(2),…,u(k)使得数据点投射到每个方向 u(i) u(i)时的误差最小。

图4 PCA实例，将2维数据点压缩为1维数据点，找到新的方向 u1 u1，使得投射误差（图中的垂线距离如 xi xi到 x˜i x~i）最小

注意：PCA和线性回归的区别，PCA是保证投射的误差（图5右的黄线）最小，而线性回归是保证沿 y y方向的误差（图5左的黄线）最小.

图5 线性回归和PCA优化目标的区别

PCA算法步骤

1. 数据预处理：mean normalization： μj=1m∑i=1mx(i)j,x(i)j=xj−μj μj=1m∑i=1mxj(i),xj(i)=xj−μj；feature scaling：(可选，不同特征范围差距过大时需要) ，  x(i)j=x(i)−μjσj xj(i)=x(i)−μjσj 

2. 计算协方差矩阵(Convariance Matrix)

Σ=1m∑i=1mx(i)(x(i))TorΣ=1mXTX Σ=1m∑i=1mx(i)(x(i))TorΣ=1mXTX

3. 计算协方差矩阵 Σ Σ的特征向量  [U, S, V] = svd(Sigma) 

4. 选择U矩阵的前k个列向量作为k个主元方向，形成矩阵 Ureduce Ureduce

5. 对于每个原始数据点 x x( x∈Rn x∈Rn)，其降维后的数据点 z z( z∈Rk z∈Rk)为  z=UTreducex z=UreduceTx

应用PCA

重构数据：对于降维后k维数据点z，将其恢复n维后的近似点为  xapporx(≈x)=Ureducez xapporx(≈x)=Ureducez

选择k值

• 平均投射误差(Average square projection error)： 1m∑i=1m||x(i)−x(i)approx||2 1m∑i=1m||x(i)−xapprox(i)||2
• total variation:  1m∑i=1m||x(i)||2 1m∑i=1m||x(i)||2
• 选择最小的k值使得  1m∑i=1m||x(i)−x(i)approx||21m∑i=1m||x(i)||2≤0.01(0.05) 1m∑i=1m||x(i)−xapprox(i)||21m∑i=1m||x(i)||2≤0.01(0.05)，也可以使用SVD分解后的S矩阵进行选择  1−∑i=1kSii∑i=1nSii≤0.01(0.05) 1−∑i=1kSii∑i=1nSii≤0.01(0.05)

应用PCA的建议

• 用于加速监督式学习：(1) 对于带标签的数据，去掉标签后进行PCA数据降维，(2)使用降维后的数据进行模型训练，(3) 对于新的数据点，先PCA降维得到降维后数据，带入模型获得预测值。注：应仅用训练集数据进行PCA降维获取映射 x(i)→z(i) x(i)→z(i)，然后将该映射（PCA选择的主元矩阵 Ureduce Ureduce）应用到验证集和测试集
• 不要用PCA阻止过拟合，用regularization。
• 在使用PCA之前，先用原始数据进行模型训练，如果不行，再考虑使用PCA；而不要上来直接使用PCA。

来看我们这次课的任务：

•数据Cat4D3Groups是4维观察数据， 
•请先采用MDS方法降维到3D，形成Cat3D3Groups数据，显示并观察。 
•对Cat3D3Groups数据采用线性PCA方法降维到2D，形成Cat2D3Groups数据，显示并观察。 

•对Cat2D3Groups数据采用K-Mean方法对数据进行分类并最终确定K，显示分类结果。 
•对Cat2D3Groups数据采用Hierarchical分类法对数据进行分类，并显示分类结果。

理论一旦推导完成，代码写起来就很轻松：

Part 1：降维处理 
MDA：

# -*- coding:gb2312 -*-from pylab import *import numpy as npfrom mpl_toolkits.mplot3d import Axes3Ddef print_D(data):
    N = np.shape(data)[0]
    d = np.zeros((N, N))    for i in range(N):
        c = data[i, :]        for j in range(N):
            e = data[j, :]
            d[i, j] = np.sqrt(np.sum(np.power(c - e, 2)))    return ddef MDS(D, K):
    N = np.shape(D)[0]
    D2 = D ** 2
    H = np.eye(N) - 1.0/N
    T = -0.5 * np.dot(np.dot(H, D2), H)
    eigVal, eigVec = np.linalg.eig(T)
    indices = np.argsort(eigVal) # 返回从小到大的索引值
    indices = indices[::-1] # 反转

    eigVal = eigVal[indices] # 特征值从大到小排列
    eigVec = eigVec[:, indices] # 排列对应特征向量

    m = eigVec[:, :K]
    n = np.diag(np.sqrt(eigVal[:K]))
    X = np.dot(m, n)    return X# test'''
data = genfromtxt("CAT4D3GROUPS.txt")
D = print_D(data)
# print D

# 4D 转 3D
CAT3D3GROUPS = MDS(D, 3)
# print CAT3D3GROUPS
# D_3D = print_D(CAT3D3GROUPS)
# print D_3D
figure(1)
ax = subplot(111,projection='3d')
ax.scatter(CAT3D3GROUPS[:, 0], CAT3D3GROUPS[:, 1], CAT3D3GROUPS[:, 2], c = 'b')
ax.set_zlabel('Z') #坐标轴
ax.set_ylabel('Y')
ax.set_xlabel('X')
title('MDS_4to3')

# 4D 转 2D
CAT2D3GROUPS = MDS(D, 2)
# print CAT2D3GROUPS
# D_2D = print_D(CAT2D3GROUPS)
# print D_2D
figure(2)
plot(CAT2D3GROUPS[:, 0], CAT2D3GROUPS[:, 1], 'b.')
xlabel('x')
ylabel('y')
title('MDS_4to2')
'''

PDA：

# -*- coding:gb2312 -*-from pylab import *from numpy import *from mpl_toolkits.mplot3d import Axes3Ddef PCA(data, K):
    # 数据标准化
    m = mean(data, axis=0) # 每列均值
    data -= m    # 协方差矩阵
    C = cov(transpose(data))    # 计算特征值特征向量，按降序排序
    evals, evecs = linalg.eig(C)
    indices = argsort(evals) # 返回从小到大的索引值
    indices = indices[::-1] # 反转

    evals = evals[indices] # 特征值从大到小排列
    evecs = evecs[:, indices] # 排列对应特征向量
    evecs_K_max = evecs[:, :K] # 取最大的前K个特征值对应的特征向量

    # 产生新的数据矩阵
    finaldata = dot(data, evecs_K_max)    return finaldata# test'''
data = genfromtxt("CAT4D3GROUPS.txt")

# 4D 转 3D
data_PCA = PCA(data, 3)
# print data_PCA
figure(1)
ax = subplot(111, projection='3d')
ax.scatter(data_PCA[:, 0], data_PCA[:, 1], data_PCA[:, 2], c='b')
ax.set_zlabel('Z') #坐标轴
ax.set_ylabel('Y')
ax.set_xlabel('X')
title('PCA_4to3')

# 4D 转 2D
data_PCA = PCA(data, 2)
print data_PCA
figure(2)
plot(data_PCA[:, 0], data_PCA[:, 1], 'b.')
xlabel('x')
ylabel('y')
title('PCA_4to2')
'''
       
       

代码里的注释啰啰嗦嗦应该解释的很清楚，这里不再赘述，看结果：

1、用MDS方法降维到3D，形成Cat3D3Groups数据： 
共两个函数，辅助函数用来生成欧氏距离矩阵，MDS函数用于降维。 


通过输出的距离矩阵可以看出，降维前后欧氏距离误差小于10^-4，证明算法有效。同时旋转3D图像也可以明显找出2D平面图的视角

2、用PCA方法降维到2D，形成Cat2D3Groups数据： 


用PCA直接对4D数据降维后的结果与MDS等价，证明算法有效。同时旋转3D图像也可以明显找出2D平面图的视角。

3、总结分析： 
先用MDS算法将4D数据降到3D，再用PCA降到2D。 

与MDS降维生成的2D图像及数据对比，误差忽略不计，证明算法有效，同时证明MDS和PCA算法在进行小批量数据降维处理上效果类似。

Part 2：聚类分析： 
数据用前面降维之后的二维数据。K-means聚类分析的程序主要参考《Machine Learning in Action》- Peter Harrington这本书第十章，我自己添加了选择最优K值的功能：

# -*- coding:gb2312 -*-import numpy as npfrom pylab import *from numpy import *# 求欧氏距离def euclDistance(vector1, vector2):
    return np.sqrt(np.sum(np.power(vector2 - vector1, 2)))# 返回某个值在列表中的全部索引值def myfind(x, y):
    return [ a for a in range(len(y)) if y[a] == x]# 初始化聚类点def initCentroids(data, k):
    numSamples, dim = data.shape
    centroids = np.zeros((k, dim))    for i in range(k):
        index = int(np.random.uniform(0, numSamples))
        centroids[i, :] = data[index, :]    return centroids# K-mean 聚类def K_mean(data, k):
    ## step 1: 初始化聚点
    centroides = initCentroids(data, k)
    numSamples = data.shape[0]
    clusterAssment = np.zeros((numSamples, 2)) # 保存每个样本点的簇索引值和误差
    clusterChanged = True

    while clusterChanged:
        clusterChanged = False
        global sum
        sum = []        # 对每一个样本点
        for i in xrange(numSamples):
            minDist = np.inf # 记录最近距离
            minIndex = 0 # 记录聚点

            ## step 2: 找到距离最近的聚点
            for j in range(k):
                distance = euclDistance(centroides[j, :], data[i, :])                if distance < minDist:
                    minDist  = distance
                    minIndex = j            ## step 3: 将该样本归到该簇
            if clusterAssment[i, 0] != minIndex:
                clusterChanged = True # 前后分类相同时停止循环
                clusterAssment[i, :] = minIndex, minDist ** 2 # 记录簇索引值和误差


        ## step 4: 更新聚点
        for j in range(k):
            index = myfind(j, clusterAssment[:, 0])
            pointsInCluster = data[index, :] # 返回属于j簇的data中非零样本的目录值，并取出样本

            centroides[j, :] = np.mean(pointsInCluster, axis=0) # 求列平均
            # 返回cost funktion值
            suml = 0
            lenght = pointsInCluster.shape[0]            for l in range(lenght):
                dis = euclDistance(centroides[j, :], pointsInCluster[l, :])
                suml += dis ** 2 / lenght
            sum.append(suml)
    cost = np.sum(sum) / k    print cost    return centroides, clusterAssment, cost# 画出分类前后结果def showCluster(data, k, centroides, clusterAssment):
    numSamples, dim = data.shape
    mark = ['r.', 'b.', 'g.', 'k.', '^r', '+r', 'sr', 'dr', '<r', 'pr']    # draw all samples
    for i in xrange(numSamples):
        markIndex = int(clusterAssment[i, 0])
        figure(2)
        plt.plot(data[i, 0], data[i, 1], mark[markIndex])
        plt.title('K-means')
    mark = ['Dr', 'Db', 'Dg', 'Dk', '^b', '+b', 'sb', 'db', '<b', 'pb']    # draw the centroids
    for i in range(k):
        figure(2)
        plt.plot(centroides[i, 0], centroides[i, 1], mark[i], markersize = 12)
       
       

其中三个辅助函数用于求欧氏距离，返回矩阵索引值和画图，k-mean函数用于聚类，当所有样本点到其所属聚类中心距离不变时，输出聚类结果，并返回cost function的值。

Cost function计算方法：对每个簇，求所有点到所属聚类中心的欧氏距离，平方后取均值E。聚类结束后，所有簇E值求和取平均得到cost function的值。

不同K值下的分类结果如下（标明聚类中心）： 

主观判断，k = 4时聚类结果最优。用Elbow方法选择K值结果如下： 

发现在K = 2时cost function值下降最为明显，与之前判断的结果不符。思考后发现，K=1时聚类没有意义，所以上图并不能有效选择K值，调整后结果如下： 

明显看出，k = 4时cost function下降极为明显。与主观判断结果相符。

Hierarchical分类，参考网上代码，出处不记得了：

# -*- coding:gb2312 -*-import numpy as npimport matplotlib.pyplot as pltimport MDSimport PCAdef yezi(clust):
    if clust.left == None and clust.right == None :        return [clust.id]    return yezi(clust.left) + yezi(clust.right)def Euclidean_distance(vector1,vector2):
    length = len(vector1)
    TSum = sum([pow((vector1[i] - vector2[i]),2) for i in range(len(vector1))])
    SSum = np.sqrt(TSum)    return SSumdef loadDataSet(fileName):
    a = []    with open(fileName, 'r') as f:
        data = f.readlines()  #txt中所有字符串读入data
        for line in data:
            odom = line.split()        #将单个数据分隔开存好
            numbers_float = map(float, odom) #转化为浮点数
            a.append(numbers_float) #print numbers_float
    a = np.array(a)    return aclass bicluster:
    def __init__(self, vec, left=None,right=None,distance=0.0,id=None):
        self.left = left
        self.right = right  #每次聚类都是一对数据，left保存其中一个数据，right保存另一个
        self.vec = vec      #保存两个数据聚类后形成新的中心
        self.id = id
        self.distance = distancedef list_array(wd, clo):
    len_=len(wd)
    xc=np.zeros([len_, clo])    for i in range(len_):
        ad = wd[i]
        xc[i, :] = ad    return xcdef hcluster(data, n) :
    [row,column] = data.shape
    data = list_array(data, column)

    biclusters = [bicluster(vec = data[i], id = i) for i in range(len(data))]
    distances = {}
    flag = None
    currentclusted = -1
    while(len(biclusters) > n) : #假设聚成n个类
        min_val = np.inf   #Python的无穷大
        biclusters_len = len(biclusters)        for i in range(biclusters_len-1) :            for j in range(i + 1, biclusters_len):                #print biclusters[i].vec

                if distances.get((biclusters[i].id,biclusters[j].id)) == None:                    #print biclusters[i].vec
                    distances[(biclusters[i].id,biclusters[j].id)] = Euclidean_distance(biclusters[i].vec,biclusters[j].vec)
                d = distances[(biclusters[i].id,biclusters[j].id)]                if d < min_val:
                    min_val = d
                    flag = (i,j)
        bic1,bic2 = flag #解包bic1 = i , bic2 = j
        newvec = [(biclusters[bic1].vec[i] + biclusters[bic2].vec[i])/2 for i in range(len(biclusters[bic1].vec))] #形成新的类中心，平均
        newbic = bicluster(newvec, left=biclusters[bic1], right=biclusters[bic2], distance=min_val, id = currentclusted) #二合一
        currentclusted -= 1
        del biclusters[bic2] #删除聚成一起的两个数据，由于这两个数据要聚成一起
        del biclusters[bic1]
        biclusters.append(newbic)#补回新聚类中心
        clusters = [yezi(biclusters[i]) for i in range(len(biclusters))] #深度优先搜索叶子节点，用于输出显示
    return biclusters,clustersdef showCluster(dataSet, k, num_mark):
    numSamples, dim = dataSet.shape
    mark = ['r.', 'b.', 'g.', 'k.', '^r', '+r', 'sr', 'dr', '<r', 'pr']    # draw all samples
    for i in xrange(numSamples):
        plt.plot(dataSet[i, 0], dataSet[i, 1], mark[num_mark])

    plt.xlabel('X')
    plt.ylabel('Y')
    plt.title('Hierarchical')if __name__ == "__main__":    # 加载数据
    dataMat =np.genfromtxt('CAT4D3GROUPS.txt')    #400*4
    dataSet = PCA.PCA(dataMat, 2)

    k,l = hcluster(dataSet, 10)  #  l返回了聚类的索引
    # 选取规模最大的k个簇，其他簇归为噪音点
    for j in range(len(l)):
        m = []        for ii in range(len(l[j])):
            m.append(l[j][ii])
        m = np.array(m)
        a = dataSet[m]
        showCluster(a,len(l),j)
    plt.show()
       
       

当一个类集合中包含多个样本点时，类与类之间的距离取Group Average：把两个集合中的点两两的欧氏距离全部放在一起求平均值，分类结果如下： 

重复运行后分类结果并未有太大变化。主观判断，从分成3类及4类的结果看，Hierarchical分类方法效果不如K-mean聚类效果好。

•