主成分分析 matlab手把手教操作、SPSS、python實例分析

更新：

　　這次決定用matlab手把手一步一步實現一遍。

令是一個的觀測矩陣，觀測向量的樣本均值M，由下式給出：

　　Matlab中我們給出矩陣.

　　　　 　　

　　用a，b度量矩陣X的行列數並分別計算每一列計算M的值

　　　　 

　　樣本均值是散列圖的中心，對，令

　　 

矩陣的列

具有零樣本均值，這樣的B稱為平均偏差形式。

　　　　 

主成分分析的目標是找到一個正交矩陣,確定一個變量代換,或並具有新的變量y₁,...,y_p兩兩無關的性質，且整理后的方差具有遞減順序。

變量的正交變換說明，每一個觀測向量X_k得到一個“新名稱”Y_k，使得，注意到Y_k是相對於以P的列為基的坐標向量。為對角形，P記為正交矩陣，它的列是單位特征向量，那么.

 

（樣本）協方差矩陣是一個矩陣S，其定義為：

　　 

由於任何具有形式的矩陣是半正定的，所以S也是半正定的。

 

數據的總方差是指S中對角線上方差的綜合。一般地，一個方陣S中對角元素之和稱為矩陣的跡，記作tr(S),這樣

{總方差}=tr{S}

S中的元素稱為x_i和x_j的協方差。 　　

　　　　

將特征值從大到小排序，取貢獻率到90%的前k個特征。

　　　　 

可得到P、D

　　　　

 下面貼完整代碼:

 

clc;
clear all;
close all;

X=[2.5 2.4;
   0.5 0.7;
   2.2 2.9;
   1.9 2.2;
   3.1 3.0;
   2.3 2.7;
   2 1.6;
   1 1.1;
   1.5 1.6;
   1.1 0.9;];
X=X'

% X=[74 87 84 88 74 86 69 73 64;
%    85 83 83 77 69 84 74 85 84;
%    83 91 89 85 87 86 83 86 85;
%    69 100 82 96 84 82 97 98 76;
%    97 48 89 36 46 53 88 89 97;
%    59 98 93 94 98 100 79 83 61;];
% X=X';

% X=[2 0 -1.4;
%    2.2 0.2 -1.5;
%    2.4 0.1 -1;
%    1.9 0 -1.2;];
% X=X';
[a,b]=size(X);
M=sum(X)/a;


for i=1:b
    B(:,i)=X(:,i)-M(i);
    %B=zscore(X);
end

S=1/(a-1)*B*B';
[vector,value]=eig(S);
vector
value=diag(value)
varine=sum(value);

[value_sort,subscript]=sort(value,'descend');
value_sort;
subscript;
value_sort=value_sort/sum(value_sort);
compare=0;
sign=0;
for i=1:b
    if compare<0.9
        sign=sign+1;
        compare=compare+value_sort(i);
    end
end

for i=1:sign
    P(:,i)=vector(:,subscript(i));
end
P

D=zeros(sign,sign);
for i=1:sign
    D(i,i)=value(subscript(i));
end
D

　　

 

---

 

　　主成分分析（Principal Component Analysis，PCA）， 是一種統計方法。通過正交變換將一組可能存在相關性的變量轉換為一組線性不相關的變量，轉換后的這組變量叫主成分。

　　主成分分析，是考察多個變量間相關性一種多元統計方法，研究如何通過少數幾個主成分來揭示多個變量間的內部結構，即從原始變量中導出少數幾個主成分，使它們盡可能多地保留原始變量的信息，且彼此間互不相關.通常數學上的處理就是將原來P個指標作線性組合，作為新的綜合指標。

　　簡單說，主成分分析的作用是降維。通過降維將原來多變量解釋的問題，映射到更少指標，轉換成少變量的可解釋性問題。但是注意經過主成分分析后的變量與原變量不存在邏輯關系，僅僅是存在線性組合的關系。[1]

。

 一、算法原理：

---

 

輸入：樣本集D={x₁,x_2...x_m}；

　　　低維空間維數d'.

過程：

1.對所有樣本進行中心化：

2.計算樣本的協方差矩陣XX^T；

3.對協方差矩陣XX^T做特征值分解；

4.取最大的d'個特征值所對應的特征向量w₁，w₂...w_d';

輸出：投影矩陣W*=（w₁，w₂...w_d'） .[2]

---

 

二、PCA原理

　　

　　

　　

 

 

 

三、SPSS進行主成分分析

　　由於SPSS本身就是一個用於數據分析的軟件，因此操作簡單無需編程，即可直觀感受主成分分析帶來的效果。

　　先胡亂編制了一些數據：

　　

　　

　　在SPSS里，點擊分析->降維->因子，在彈出的對話框中，將需要分析的變量都送入變量欄中。根據個人需要在描述、提取、旋轉、得分、選項中勾選。此處我們注意在提取中勾選主成分。

　　

 

點擊“確定”： 

　　

　　

最后我們可以看到提取了兩個主成分

　　

 

　　

觀察兩個主成分中的貢獻率，我們會發現第一個主成分包含貢獻率較高的項為數學、物理、化學、生物，實際意義即理科，第二主成分包含歷史、地理，即文科。具有良好解釋性。

 

四、python代碼實現主成分分析

pca.py

# -*- coding: utf-8 -*-
"""
Created on Sun Feb 28 10:04:26 2016
PCA source code
@author: liudiwei
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

#計算均值,要求輸入數據為numpy的矩陣格式，行表示樣本數，列表示特征    
def meanX(dataX):
    return np.mean(dataX,axis=0)#axis=0表示按照列來求均值，如果輸入list,則axis=1


#計算方差,傳入的是一個numpy的矩陣格式，行表示樣本數，列表示特征    
def variance(X):
    m, n = np.shape(X)
    mu = meanX(X)
    muAll = np.tile(mu, (m, 1))    
    X1 = X - muAll
    variance = 1./m * np.diag(X1.T * X1)
    return variance

#標准化,傳入的是一個numpy的矩陣格式，行表示樣本數，列表示特征    
def normalize(X):
    m, n = np.shape(X)
    mu = meanX(X)
    muAll = np.tile(mu, (m, 1))    
    X1 = X - muAll
    X2 = np.tile(np.diag(X.T * X), (m, 1))
    XNorm = X1/X2
    return XNorm

"""
參數：
	- XMat：傳入的是一個numpy的矩陣格式，行表示樣本數，列表示特征    
	- k：表示取前k個特征值對應的特征向量
返回值：
	- finalData：參數一指的是返回的低維矩陣，對應於輸入參數二
	- reconData：參數二對應的是移動坐標軸后的矩陣
"""  
def pca(XMat, k):
    average = meanX(XMat) 
    m, n = np.shape(XMat)
    data_adjust = []
    avgs = np.tile(average, (m, 1))
    data_adjust = XMat - avgs
    covX = np.cov(data_adjust.T)   #計算協方差矩陣
    featValue, featVec=  np.linalg.eig(covX)  #求解協方差矩陣的特征值和特征向量
    index = np.argsort(-featValue) #按照featValue進行從大到小排序
    finalData = []
    if k > n:
        print("k must lower than feature number")
        return
    else:
        #注意特征向量時列向量，而numpy的二維矩陣(數組)a[m][n]中，a[1]表示第1行值
        selectVec = np.matrix(featVec.T[index[:k]]) #所以這里需要進行轉置
        finalData = data_adjust * selectVec.T 
        reconData = (finalData * selectVec) + average  
    return finalData, reconData

def loaddata(datafile):
    return np.array(pd.read_csv(datafile,sep="\t",header=-1)).astype(np.float)


def plotBestFit(data1, data2):	  
    dataArr1 = np.array(data1)
    dataArr2 = np.array(data2)
    
    m = np.shape(dataArr1)[0]
    axis_x1 = []
    axis_y1 = []
    axis_x2 = []
    axis_y2 = []
    for i in range(m):
        axis_x1.append(dataArr1[i,0])
        axis_y1.append(dataArr1[i,1])
        axis_x2.append(dataArr2[i,0]) 
        axis_y2.append(dataArr2[i,1])				  
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(axis_x1, axis_y1, s=50, c='red', marker='s')
    ax.scatter(axis_x2, axis_y2, s=50, c='blue')
    plt.xlabel('x1'); plt.ylabel('x2');
    plt.savefig("outfile.png")
    plt.show()	

#簡單測試
#數據來源：http://www.cnblogs.com/jerrylead/archive/2011/04/18/2020209.html
def test():
    X = [[2.5, 0.5, 2.2, 1.9, 3.1, 2.3, 2, 1, 1.5, 1.1],
         [2.4, 0.7, 2.9, 2.2, 3.0, 2.7, 1.6, 1.1, 1.6, 0.9]]
    XMat = np.matrix(X).T  
    k = 2
    return pca(XMat, k)

#根據數據集data.txt
def main():    
    datafile = "data.txt"
    XMat = loaddata(datafile)
    k = 2
    return pca(XMat, k)
    
if __name__ == "__main__":
    finalData, reconMat = main()
    plotBestFit(finalData, reconMat)

 

data.txt

5.1	3.5	1.4	0.2
4.9	3.0	1.4	0.2
4.7	3.2	1.3	0.2
4.6	3.1	1.5	0.2
5.0	3.6	1.4	0.2
5.4	3.9	1.7	0.4
4.6	3.4	1.4	0.3
5.0	3.4	1.5	0.2
4.4	2.9	1.4	0.2
4.9	3.1	1.5	0.1
5.4	3.7	1.5	0.2
4.8	3.4	1.6	0.2
4.8	3.0	1.4	0.1
4.3	3.0	1.1	0.1
5.8	4.0	1.2	0.2
5.7	4.4	1.5	0.4
5.4	3.9	1.3	0.4
5.1	3.5	1.4	0.3
5.7	3.8	1.7	0.3
5.1	3.8	1.5	0.3
5.4	3.4	1.7	0.2
5.1	3.7	1.5	0.4
4.6	3.6	1.0	0.2
5.1	3.3	1.7	0.5
4.8	3.4	1.9	0.2
5.0	3.0	1.6	0.2
5.0	3.4	1.6	0.4
5.2	3.5	1.5	0.2
5.2	3.4	1.4	0.2
4.7	3.2	1.6	0.2
4.8	3.1	1.6	0.2
5.4	3.4	1.5	0.4
5.2	4.1	1.5	0.1
5.5	4.2	1.4	0.2
4.9	3.1	1.5	0.1
5.0	3.2	1.2	0.2
5.5	3.5	1.3	0.2
4.9	3.1	1.5	0.1
4.4	3.0	1.3	0.2
5.1	3.4	1.5	0.2
5.0	3.5	1.3	0.3
4.5	2.3	1.3	0.3
4.4	3.2	1.3	0.2
5.0	3.5	1.6	0.6
5.1	3.8	1.9	0.4
4.8	3.0	1.4	0.3
5.1	3.8	1.6	0.2
4.6	3.2	1.4	0.2
5.3	3.7	1.5	0.2
5.0	3.3	1.4	0.2
7.0	3.2	4.7	1.4
6.4	3.2	4.5	1.5
6.9	3.1	4.9	1.5
5.5	2.3	4.0	1.3
6.5	2.8	4.6	1.5
5.7	2.8	4.5	1.3
6.3	3.3	4.7	1.6
4.9	2.4	3.3	1.0
6.6	2.9	4.6	1.3
5.2	2.7	3.9	1.4
5.0	2.0	3.5	1.0
5.9	3.0	4.2	1.5
6.0	2.2	4.0	1.0
6.1	2.9	4.7	1.4
5.6	2.9	3.6	1.3
6.7	3.1	4.4	1.4
5.6	3.0	4.5	1.5
5.8	2.7	4.1	1.0
6.2	2.2	4.5	1.5
5.6	2.5	3.9	1.1
5.9	3.2	4.8	1.8
6.1	2.8	4.0	1.3
6.3	2.5	4.9	1.5
6.1	2.8	4.7	1.2
6.4	2.9	4.3	1.3
6.6	3.0	4.4	1.4
6.8	2.8	4.8	1.4
6.7	3.0	5.0	1.7
6.0	2.9	4.5	1.5
5.7	2.6	3.5	1.0
5.5	2.4	3.8	1.1
5.5	2.4	3.7	1.0
5.8	2.7	3.9	1.2
6.0	2.7	5.1	1.6
5.4	3.0	4.5	1.5
6.0	3.4	4.5	1.6
6.7	3.1	4.7	1.5
6.3	2.3	4.4	1.3
5.6	3.0	4.1	1.3
5.5	2.5	4.0	1.3
5.5	2.6	4.4	1.2
6.1	3.0	4.6	1.4
5.8	2.6	4.0	1.2
5.0	2.3	3.3	1.0
5.6	2.7	4.2	1.3
5.7	3.0	4.2	1.2
5.7	2.9	4.2	1.3
6.2	2.9	4.3	1.3
5.1	2.5	3.0	1.1
5.7	2.8	4.1	1.3
6.3	3.3	6.0	2.5
5.8	2.7	5.1	1.9
7.1	3.0	5.9	2.1
6.3	2.9	5.6	1.8
6.5	3.0	5.8	2.2
7.6	3.0	6.6	2.1
4.9	2.5	4.5	1.7
7.3	2.9	6.3	1.8
6.7	2.5	5.8	1.8
7.2	3.6	6.1	2.5
6.5	3.2	5.1	2.0
6.4	2.7	5.3	1.9
6.8	3.0	5.5	2.1
5.7	2.5	5.0	2.0
5.8	2.8	5.1	2.4
6.4	3.2	5.3	2.3
6.5	3.0	5.5	1.8
7.7	3.8	6.7	2.2
7.7	2.6	6.9	2.3
6.0	2.2	5.0	1.5
6.9	3.2	5.7	2.3
5.6	2.8	4.9	2.0
7.7	2.8	6.7	2.0
6.3	2.7	4.9	1.8
6.7	3.3	5.7	2.1
7.2	3.2	6.0	1.8
6.2	2.8	4.8	1.8
6.1	3.0	4.9	1.8
6.4	2.8	5.6	2.1
7.2	3.0	5.8	1.6
7.4	2.8	6.1	1.9
7.9	3.8	6.4	2.0
6.4	2.8	5.6	2.2
6.3	2.8	5.1	1.5
6.1	2.6	5.6	1.4
7.7	3.0	6.1	2.3
6.3	3.4	5.6	2.4
6.4	3.1	5.5	1.8
6.0	3.0	4.8	1.8
6.9	3.1	5.4	2.1
6.7	3.1	5.6	2.4
6.9	3.1	5.1	2.3
5.8	2.7	5.1	1.9
6.8	3.2	5.9	2.3
6.7	3.3	5.7	2.5
6.7	3.0	5.2	2.3
6.3	2.5	5.0	1.9
6.5	3.0	5.2	2.0
6.2	3.4	5.4	2.3
5.9	3.0	5.1	1.8

　　代碼運行結果：

 

參考文獻：

[1] https://baike.baidu.com/item/%E4%B8%BB%E6%88%90%E5%88%86%E5%88%86%E6%9E%90/829840?fr=aladdin

[2]周志華，機器學習，清華大學出版社，2016年1月1版.

 

2019-03-09

00:05:40