主成分分析,即Principal Component Analysis(PCA),是多元统计中的重要内容,也广泛应用于机器学习和其它领域。它的主要作用是对高维数据进行降维。PCA把原先的n个特征用数目更少的k个特征取代,新特征是旧特征的线性组合,这些线性组合最大化样本方差,尽量使新的k个特征互不相关。关于PCA的更多介绍,请参考:https://en.wikipedia.org/wiki/Principal_component_analysis.

PCA的主要算法如下:

组织数据形式,以便于模型使用;

计算样本每个特征的平均值;

每个样本数据减去该特征的平均值(归一化处理);

求协方差矩阵;

找到协方差矩阵的特征值和特征向量;

对特征值和特征向量重新排列(特征值从大到小排列);

对特征值求取累计贡献率;

对累计贡献率按照某个特定比例选取特征向量集的子集合;

对原始数据(第三步后)进行转换。

其中协方差矩阵的分解可以通过按对称矩阵的特征向量来,也可以通过分解矩阵的SVD来实现,而在Scikit-learn中,也是采用SVD来实现PCA算法的。关于SVD的介绍及其原理,可以参考:矩阵的奇异值分解(SVD)(理论)。

本文将用三种方法来实现PCA算法,一种是原始算法,即上面所描述的算法过程,具体的计算方法和过程,可以参考:A tutorial on Principal Components Analysis, Lindsay I Smith. 一种是带SVD的原始算法,在Python的Numpy模块中已经实现了SVD算法,并且将特征值从大从小排列,省去了对特征值和特征向量重新排列这一步。最后一种方法是用Python的Scikit-learn模块实现的PCA类直接进行计算,来验证前面两种方法的正确性。

用以上三种方法来实现PCA的完整的Python如下:

import numpy as np
from sklearn.decomposition import PCA
import sys
#returns choosing how many main factors
def index_lst(lst, component=0, rate=0):
#component: numbers of main factors
#rate: rate of sum(main factors)/sum(all factors)
#rate range suggest: (0.8,1)
#if you choose rate parameter, return index = 0 or less than len(lst)
if component and rate:
print('Component and rate must choose only one!')
sys.exit(0)
if not component and not rate:
print('Invalid parameter for numbers of components!')
sys.exit(0)
elif component:
print('Choosing by component, components are %s......'%component)
return component
else:
print('Choosing by rate, rate is %s ......'%rate)
for i in range(1, len(lst)):
if sum(lst[:i])/sum(lst) >= rate:
return i
return 0
def main():
# test data
mat = [[-1,-1,0,2,1],[2,0,0,-1,-1],[2,0,1,1,0]]
# simple transform of test data
Mat = np.array(mat, dtype='float64')
print('Before PCA transforMation, data is:\n', Mat)
print('\nMethod 1: PCA by original algorithm:')
p,n = np.shape(Mat) # shape of Mat
t = np.mean(Mat, 0) # mean of each column
# substract the mean of each column
for i in range(p):
for j in range(n):
Mat[i,j] = float(Mat[i,j]-t[j])
# covariance Matrix
cov_Mat = np.dot(Mat.T, Mat)/(p-1)
# PCA by original algorithm
# eigvalues and eigenvectors of covariance Matrix with eigvalues descending
U,V = np.linalg.eigh(cov_Mat)
# Rearrange the eigenvectors and eigenvalues
U = U[::-1]
for i in range(n):
V[i,:] = V[i,:][::-1]
# choose eigenvalue by component or rate, not both of them euqal to 0
Index = index_lst(U, component=2) # choose how many main factors
if Index:
v = V[:,:Index] # subset of Unitary matrix
else: # improper rate choice may return Index=0
print('Invalid rate choice.\nPlease adjust the rate.')
print('Rate distribute follows:')
print([sum(U[:i])/sum(U) for i in range(1, len(U)+1)])
sys.exit(0)
# data transformation
T1 = np.dot(Mat, v)
# print the transformed data
print('We choose %d main factors.'%Index)
print('After PCA transformation, data becomes:\n',T1)
# PCA by original algorithm using SVD
print('\nMethod 2: PCA by original algorithm using SVD:')
# u: Unitary matrix, eigenvectors in columns
# d: list of the singular values, sorted in descending order
u,d,v = np.linalg.svd(cov_Mat)
Index = index_lst(d, rate=0.95) # choose how many main factors
T2 = np.dot(Mat, u[:,:Index]) # transformed data
print('We choose %d main factors.'%Index)
print('After PCA transformation, data becomes:\n',T2)
# PCA by Scikit-learn
pca = PCA(n_components=2) # n_components can be integer or float in (0,1)
pca.fit(mat) # fit the model
print('\nMethod 3: PCA by Scikit-learn:')
print('After PCA transformation, data becomes:')
print(pca.fit_transform(mat)) # transformed data
main()

运行以上代码,输出结果为:


这说明用以上三种方法来实现PCA都是可行的。这样我们就能理解PCA的具体实现过程啦~~有兴趣的读者可以用其它语言实现一下哈~~

参考文献:

PCA 维基百科: https://en.wikipedia.org/wiki/Principal_component_analysis.

讲解详细又全面的PCA教程: A tutorial on Principal Components Analysis, Lindsay I Smith.


Scikit-learn的PCA介绍:http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html.

Python使用三种方法实现PCA算法[转]

主成分分析(PCA) vs 多元判别式分析(MDA) PCA和MDA都是线性变换的方法,二者关系密切.在PCA中,我们寻找数据集中最大化方差的成分,在MDA中,我们对类间最大散布的方向更感兴趣. 一句 ...

python字符串连接的三种方法及其效率、适用场景详解

python字符串连接的方法,一般有以下三种:方法1:直接通过加号(+)操作符连接website=& 39;python& 39;+& 39;tab& 39;+&amp ...

python每次处理一个字符的三种方法

python每次处理一个字符的三种方法 a_string = "abccdea" print 'the first' for c in a_string: print ord(c) ...

python更新数据库脚本三种方法

最近项目的两次版本迭代中,根据业务需求的变化,需要对数据库进行更新,两次分别使用了不同的方式进行更新. 第一种:使用python的MySQLdb模块利用原生的sql语句进行更新 import MySQ ...

python下载文件的三种方法

Python开发中时长遇到要下载文件的情况,最常用的方法就是通过Http利用urllib或者urllib2模块. 当然你也可以利用ftplib从ftp站点下载文件.此外Python还提供了另外一种方法 ...

python网络编程调用recv函数完整接收数据的三种方法

最近在使用python进行网络编程开发一个通用的tcpclient测试小工具.在使用socket进行网络编程中,如何判定对端发送一条报文是否接收完成,是进行socket网络开发必须要考虑的一个问题.这 ...

python 多线程编程之threading模块(Thread类)创建线程的三种方法

摘录 python核心编程 上节介绍的thread模块,是不支持守护线程的.当主线程退出的时候,所有的子线程都将终止,不管他们是否仍在工作. 本节开始,我们开始介绍python的另外多线程模块thre ...

python—字符串拼接三种方法

python—字符串拼接三种方法   1.使用加号(+)号进行拼接 字符串拼接直接进行相加就可以,比较容易理解,但是一定要记得,变量直接相加,不是变量就要用引号引起来,不然会出错,另外数字是要转换为字 ...

随机推荐

IOS数据存储之Sqlite数据库

前言: 之前学习了数据存储的NSUserDefaults,归档和解档,沙盒文件存储,但是对于数据量比较大,需要频繁查询,删除,更新等操作的时候无论从效率上还是性能上,上述三种明显不能满足我们的日常开发 ...

架构实例之Demo_JSP_JavaBean

架构实例之Demo_JSP_JavaBean 1.开发工具和开发环境      开发工具: MyEclipse10,JDK1.6.0_13(32位),Tomcat7.0(32位),mysql5.7.1 ...

mongo-c-driver使用VS2013编译

1.下载mongo-c-driver源码文件 使用github来下载. git clone https://github.com/mongodb/mongo-c-driver.git 下载完之后,进入 ...

MIT 6.828 JOS学习笔记1. Lab 1 Part 1: PC Bootstrap

Lab 1: Booting a PC Part 1: PC Bootstrap 介绍这一部分知识的目的就是让你能够更加熟悉x86汇编语言,以及PC启动的整个过程,而且也会首次学习使用QEMU软件来仿 ...

Traffic Lights

Traffic Lights time limit per test 2 seconds memory limit per test 256 megabytes input standard inpu ...

nginx location 匹配顺序

location 匹配的原型是这样的:location [=|~|~*|^~|@] /uri/ { … } “=”是精确匹配“@”是命名的location ,在正常的location 匹配中不会使用, ...

URL重写以后发布到IIS找不到页面

1.读取必须勾选,否则无法加载资源文件(img,css等) c:\windows\microsoft.net\framework\v2.0.50727\aspnet_isapi.dll

UVa 1467 (贪心+暴力) Installations

题意: 一共有n项服务,每项服务有安装的时间s和截止时间d.对于每项任务,如果有一项超出截止时间,惩罚值为所超出时间的长度.问如何安装才能使惩罚值最大的两个任务的惩罚值之和最小. 分析: 如果是求总惩 ...

【2011 Greater New York Regional 】Problem B	The Rascal Triangle

一个简单的规律题,每一列都是一个等差数列: 代码: #include #define ll long long using namespace std; int main( ...