參考《機器學習實戰》，使用MNIST資料集

首先給出PCA的演算法

輸入：樣本集 

           低維空間維數 d‘

過程：

1. 對所有樣本進行中心化：
2. 計算樣本的協方差矩陣
3. 對協方差矩陣做特徵值分解
4. 取最大的d'個特徵值所對應的特徵向量

輸出：投影矩陣

# -*- coding: utf-8 -*-
"""
Created on Sun Jul 22 15:38:03 2018

@author: Li Qingquan
"""

import chainer
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D

# with reference to "Machine Learning in Action"
def PCA(dataMat, topNfeat=9999999):
    meanVals = np.mean(dataMat, axis = 0)  # calculate mean
    meanRemoved = dataMat - meanVals  # minus mean
    covMat = np.cov(meanRemoved, rowvar=0)  # calculate the covariance matrix
    eigVals, eigVects = np.linalg.eig(np.mat(covMat))  # calculate eigenvalues
    # sort eigenvalues
    eigValInd = np.argsort(eigVals)
    eigValInd = eigValInd[: -(topNfeat + 1) : -1]
    redEigVects = eigVects[:, eigValInd]
    # dimensionality reduction
    lowDDataMat = meanRemoved * redEigVects
    reconMat = (lowDDataMat * redEigVects.T) + meanVals
    return lowDDataMat, reconMat

if __name__ == "__main__":
    train, test = chainer.datasets.get_mnist()
    img = []
    label = []
    for i in range(10000):  # get numbers 3, 6, 8
        if train[i][1] == 3 or train[i][1] == 6 or train[i][1] == 8:
            img.append(train[i][0])
            label.append(train[i][1])
    lowDDataMat, reconMat = PCA(np.array(img), 3)
    # draw
    fig = plt.figure()
    axes = fig.add_subplot(111, projection='3d')
    for i in range(len(label)):
        if label[i] == 3:
            axes.scatter(np.array(lowDDataMat[i, 0]), np.array(lowDDataMat[i, 1]), np.array(lowDDataMat[i, 2]), color='red')
        if label[i] == 6:
            axes.scatter(np.array(lowDDataMat[i, 0]), np.array(lowDDataMat[i, 1]), np.array(lowDDataMat[i, 2]), color='green')
        if label[i] == 8:
            axes.scatter(np.array(lowDDataMat[i, 0]), np.array(lowDDataMat[i, 1]), np.array(lowDDataMat[i, 2]), color='blue')
    plt.title("PCA")
    plt.show()