python降维方法_Python实现PCA降维算法及其应用

本帖最后由 Oner 于 2017-8-10 13:26 编辑

问题导读：

1. PCA 降维算法原理是什么？

2. 如何使用Python实现PCA降维算法？

3. PCA降维后的数据常用于什么场景？

4. Kmeans算法原理是什么？

5. 如何使用Spark Mllib自带的Kmeans算法?

前言首先我们先确定一个事实，就是我们在做ML(机器学习)的时候，绝不是算法第一的。我们在很多时候选择一个或者说决定一个模型开始训练，我们首先应该想的是：

数据来源(ETL的T)

数据的格式化(ETL的T)

数据采用的训练模型

模型的展示PCA降维算法

PCA简单的说，它是一种通用的降维工具。在我们处理高维数据的时候，

可以(原文“了能”)降低后续计算的复杂度，在“预处理”阶段通常要先对原始数据进行降维，

而PCA就是干这个事的 本质上讲，PCA就是将高维的数据通过线性变换投影到低维空间上去

具体的数学原理我推荐这个(http://blog.csdn.net/xiaojidan2011/article/details/11595869)

降维当然意味着信息的丢失，不过鉴于实际数据本身常常存在的相关性，我们可以想办法在降维的同时将信息的损失尽量降低。总结一下PCA的算法步骤：

设有m条n维数据。

1)将原始数据按列组成n行m列矩阵X

2)将X的每一行(代表一个属性字段)进行零均值化，即减去这一行的均值

3)求出协方差矩阵C=

QQ截图20170810101305.png (682 Bytes, 下载次数: 2)

2017-8-10 10:13 上传

4)求出协方差矩阵的特征值及对应的特征向量

5)将特征向量按对应特征值大小从上到下按行排列成矩阵，取前k行组成矩阵P

6)Y=PX即为降维到k维后的数据

我们用python试着实现下：

[mw_shl_code=python,true]# -*- coding: utf-8 -*-

from math import *

import random as rd

import numpy as np

import matplotlib as mpl

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

def zeroMean(dataMat):

meanVal = np.mean(dataMat,axis = 0)#计算该轴上的统计值(0为列，1为行)

newData = dataMat - meanVal

return newData,meanVal

def pca(dataMat,percent=0.99):

'''求协方差矩阵

若rowvar=0，说明传入的数据一行代表一个样本，若非0

说明传入的数据一列代表一个样本。因为newData每一行代表一个样本，所以将rowvar设置为0 '''

newData,meanVal=zeroMean(dataMat)

covMat=np.cov(newData,rowvar=0)

eigVals,eigVects = np.linalg.eig(np.mat(covMat))

n=percentage2n(eigVals,percent)          #要达到percent的方差百分比，需要前n个特征向量

print str(n) + u"vectors"

eigValIndice=np.argsort(eigVals)            #对特征值从小到大排序

n_eigValIndice=eigValIndice[-1:-(n+1):-1]   #最大的n个特征值的下标

n_eigVect=eigVects[:,n_eigValIndice]        #最大的n个特征值对应的特征向量

lowDDataMat=newData * n_eigVect               #低维特征空间的数据

reconMat=(lowDDataMat * n_eigVect.T) + meanVal  #重构数据

return reconMat,lowDDataMat,n

def percentage2n(eigVals,percentage):

sortArray=np.sort(eigVals)   #升序

sortArray=sortArray[-1::-1]  #逆转，即降序

arraySum=sum(sortArray)

tmpSum=0

num=0

for i in sortArray:

tmpSum += i

num += 1

if tmpSum >= arraySum * percentage:

return num

if __name__ == '__main__':

data = np.random.randint(1,10,size = (3,5))

fig = plt.figure()

ax = plt.subplot(111,projection='3d')

#ax.scatter(data[0],data[1],data[2],c='y') #绘制数据点

ax.set_zlabel('Z') #坐标轴

ax.set_ylabel('Y')

ax.set_xlabel('X')

#plt.show()

print data

fin = pca(data,0.9)

mat =fin[1]

print mat

#ax.scatter(mat[0],mat[1],mat[2],c='y') #绘制数据点

#plt.show()[/mw_shl_code]

应用

我们对PCA后降维后的数据最直接的应用是聚类，这里我们还是选择kmeans算法：

kmeans

具体算法可百度，我这里提下不同类型变量相异度计算方法：

标量也就是无方向意义的数字，也叫标度变量:

一种很自然的想法是用两者的欧几里得距离来作为相异度，欧几里得距离的定义如下：

QQ截图20170810104218.png (2.17 KB, 下载次数: 3)

2017-8-10 10:42 上传

对于向量，由于它不仅有大小而且有方向，所以闵可夫斯基距离不是度量其相异度的好办法，一种流行的做法是用两个向量的余弦度量，其度量公式为：

QQ截图20170810104336.png (1.43 KB, 下载次数: 1)

2017-8-10 10:43 上传k均值算法的计算过程非常直观：

1、从D中随机取k个元素，作为k个簇的各自的中心。

2、分别计算剩下的元素到k个簇中心的相异度，将这些元素分别划归到相异度最低的簇。

3、根据聚类结果，重新计算k个簇各自的中心，计算方法是取簇中所有元素各自维度的算术平均数。

4、将D中全部元素按照新的中心重新聚类。

5、重复第4步，直到聚类结果不再变化。

6、将结果输出。

我们上代码，是基于spark MLlib的：[mw_shl_code=scala,true]import org.apache.spark.SparkConf

import org.apache.spark.SparkContext

import org.apache.spark.mllib.clustering.KMeans

import org.apache.spark.mllib.linalg.Matrix

import org.apache.spark.mllib.linalg.Vectors

import org.apache.spark.mllib.linalg.distributed.RowMatrix

import scala.util.Random

object PCA {

def getRandom(num: Int) = {

(for(i

}

def main(args: Array[String]): Unit = {

val line = 5

val conf =new SparkConf().setAppName("PCA");

val sc = new SparkContext(conf)

val data =(for(i

//我们产生一个5X10的矩阵

val dataRDD = sc.parallelize(data, 2)

//RowMatrix 分布式矩阵,RowMatrix.numRows，RowMatrix.numCols

val mat: RowMatrix = new RowMatrix(dataRDD)

val pc: Matrix = mat.computePrincipalComponents(3)//PCA只需要保留前3个特征

//得到的矩阵结果

val projected: RowMatrix = mat.multiply(pc)

val newdateRDD = projected.rows

// 打印出降维的数据

val numIterations = 20//迭代的次数

val ks:Array[Int] = Array(2,3,4,5)

ks.foreach(cluster => {

val model = KMeans.train(newdateRDD, cluster,numIterations,1)

val ssd = model.computeCost(newdateRDD)

println(" when k=" + cluster + " -> "+ ssd)

})

val Knum = 3  //将目标数据分成几类

//将参数，和训练数据传入，形成模型

val clusters1 = KMeans.train(dataRDD, Knum , numIterations)//训练原始数据

val clusters2 = KMeans.train(newdateRDD, Knum , numIterations)//训练降维数据

val collect = projected.rows.collect()

println("主成分投影列矩阵:")

collect.foreach { vector => println(vector) }

val source = mat.rows.collect()

println("原始投影列矩阵:")

source.foreach { vector => println(vector) }

/*println("原始数据中心分布:")

for (c

println("  " + c.toString)

}

println("训练降维数据中心分布:")

for (c

println("  " + c.toString)

}*/

for( i

println("训练原始数据:" + source(i).toString + "属于" + clusters1.predict(source(i)).toString + "类")

}

for( i

println("训练降维数据:" + collect(i).toString + "属于" + clusters2.predict(collect(i)).toString + "类")

}

}

}[/mw_shl_code]

后记

大概就这样，最近遇到一个问题，就是数据的维度太多，特征向量都不知道选什么。很是头疼~

作者：michaelgbw

來源：http://www.jianshu.com/p/16f835a94444