盒子
盒子

主成分分析的两种语言实现

2015.10.21 nichunquan machine learning  热度 ℃

###1.主成分分析计算步骤

实现过程如下：

（1）计算标准化的值z-score（z=(x-mu)/delta）关于为什么要用z-score做标准化的原因是：减去均值等同于坐标的移动，把原始数据点的重心移到和原点重合，这样利于很多表达，比如数据的协方差矩阵可以写成XX’，若没有减去均值，则XX‘后面还要减去一些东西（还不明白可以参考多元统计分析的书）。除以标准差是为了统一并消除量纲。一个矩阵中有多个向量，有些可能表示了长度，有些表示了重量，除以标准差，才能让它们仅以“数”的概念一起比较运算。

（2）计算z-score的协方差矩阵，根据协方差矩阵计算特征值和特征向量

（3）将特征值进行有大到小排列，保留所有特征值中的前N个特征值，并保留其对应的特征向量

（4）将获得的N个特征向量映射到新的空间

###2. python

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
import numpy as np
from scipy import linalg
from sklearn.datasets import load_iris
from numpy import * 
def qr(A):
    n, m = A.shape
    R = A.copy()
    Q = np.eye(n)
    for k in range(m-1):
        x = np.zeros((n, 1))
        x[k:, 0] = R[k:, k]
        v = x
        v[k] = x[k] + np.sign(x[k,0]) * np.linalg.norm(x)
        s = np.linalg.norm(v)
        if s != 0:
            u = v / s
            R -= 2 * np.dot(u, np.dot(u.T, R))
            Q -= 2 * np.dot(u, np.dot(u.T, Q))
    Q = Q.T
    return Q, R
def eigen_qr(A):
    T = 1000
    A_copy = A.copy()
    r, c = A_copy.shape
    V = np.random.random_sample((r, r))
    for i in range(T):
        Q, _ = qr(V)
        V = np.dot(A_copy, Q)
    Q, R = qr(V)
    return R.diagonal(), Q
def my_pca(A, topNfeat=5):
    meanVals = mean(A, axis=0)
    meanRemoved = A - meanVals
    stded = meanRemoved / std(A,axis=0)
    covA=cov(stded, rowvar=0)
    eigVals, eigVects=eigen_qr(covA)
    eigValInd = argsort(eigVals)
    eigValInd = eigValInd[:-(topNfeat + 1):-1]
    Lambda = np.array(eigVals)[eigValInd]**0.5
    Q = eigVects[:, eigValInd]
    return   Lambda, Q
iris = load_iris()
my_pca(iris.data)

国内很多网站有类似的PCA用python实现的代码，但里面有些小错误，std()这个函数如果axis不加，那么它默认对整个矩阵A求标准差，结果是一个值，如果是加上axis=0那么就是对列求标准差。

###3. R语言

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
myqr <- function(A)
{
    n = nrow(A)
    m = ncol(A)
    R = A
    Q = diag(n)
    
    for (k in 1:(m-1))
    {
        x = matrix(rep(0, n), nrow = n)
        x[k:n, 1] = R[k:n, k]
        g = sqrt(sum(x^2))
        v = x
        v[k] = x[k] + sign(x[k,1])*g
        
        s = sqrt(sum(v^2))
        if (s != 0)
        {
            u = v / s
            R = R - 2 * u %*% t(u) %*% R
            Q = Q - 2 * u %*% t(u) %*% Q
        }
    }
    result <- list(t(Q), R)
    names(result) <- c("Q", "R")
    result
}
myeigen <- function(A)     #求特征值和特征向量
{
    T = 1000
    n = nrow(A)
    V = matrix(rnorm(n*n), nrow=n)
    
    for (i in 1:T)
    {
        V = myqr(V)$Q
        V = A %*% V
    }
    
    B = myqr(V)
    
    result <- list(B$Q, diag(B$R))
    names(result) <- c("vectors", "values")
    result
}
my_pca <- function(A,topNfeat)
{
  sd=scale(A,center=TRUE,scale=TRUE)
  covA=cov(sd)
  eigen=myeigen(covA)
  result <- list(t(eigen$vector[1:topNfeat,]),eigen$values[1:topNfeat]^0.5) #比较标准差，模仿自带的prcomp函数
  names(result) <- c("Q","lambda")
  result
}
iris_data_frame = data.frame(iris[,1:4])
kk=my_pca(iris_data_frame,4)

两个代码思路是一样的一些矩阵操作依据语言不同而略显差异。