主成份分析（Principal Component Analysis, PCA）最主要是利用數學的方法，將複雜的事情簡化之。主成份分析在 100 年前由英國數學家卡爾·皮爾森發明，是一個至今仍在機器學習與統計學領域中被廣泛用來分析資料、降低數據維度以及去關聯的線性降維方法。

說穿了，主成份分析就是要把二維平面降成一維的線。
而三維立體空間呢，則降成二維平面，甚至也可以再降成一條線。
再進階一點，把四維降成三維，二維，或一維。100維也可以降成一維，二維，三維，隨便你高興的降。

那麼，怎麼降呢?，應該有一定的SOP吧!! 對，這個SOP就是最短距離，也就是投影距離。

OK, 主成分分析就是去除不想要知道的成份，只留下想知道的資訊。本篇只說明PCA的基本原理，並不去討論它的演算法，因為演算法需要有高等數學的涵養才看得懂.

此篇雖不說明其演算法，但相關基礎數學還是必需知道，請先研讀基礎數學這篇說明

模型應用

假設我們要研究所有鳥類，但我們的研究的需求不同，就有不同的解法(模型)

1. 想知道鳥類分佈區域 :

此時就使用 PCA 降維，將所有的鳥類都打死，然後看地面的分佈狀況。這樣就會知道老鷹最多是在非洲，麻雀大多在台灣。(因為我不想知道老鷹會飛多高，麻雀有多膽小)

2. 想知道有幾種鳥類 :

使用 SVM 模型，叫所有的鳥類飛到個自的高度，然後依高度分類。50公尺以下是麻雀，1000公尺以上是老鷹。

3. k-mean

上述 SVM 模型，最後還是需要使用 k-mean 進行分類作業，這也是大多數人不知道的一件事。

產生資料

首先產生20組數字資料以便方便日後說明，底下為產生20組數字的程式.

import numpy as np
import pylab as plt
# 只顯示小數點下兩位數字，但實際還是以原始值計算
np.set_printoptions(precision=2)

#偽隨機產生器，每次重頭開始執行其值都一樣
rng = np.random.RandomState(1)
W = rng.rand(2, 2)

#產生二列20行的隨機數, scale愈大愈矮胖，愈小愈高瘦
X_normal = rng.normal(scale=5, size=(2, 20))

#將W @ X_normal 是二矩陣相乘，[[1,2,3],[2,2,2]] @ [[1,2],[3,4],[5,6]] = [[22,28],[18,24]]
X_orig = W @ X_normal

#計算每列的平均值，再轉成二維陣列
X_mean = X_orig.mean(axis=1).reshape(2,1)

X = X_orig - X_mean
print(X)

plt.xlim(-8,8)
plt.ylim(-8,8)
plt.scatter(X[0], X[1])
plt.show()

結果:
[[ 2.89 0.32 5.8 -6.52 3.94 -4.21 0.45 2.14 1.3 -4.98 -2.4 -3.1
0.69 -1.59 -3.64 -0.24 6.81 4.63 -2.24 -0.06]
[ 1.52 0.91 1.52 -0.88 -0.03 -1.26 -0.25 0.96 -0.89 -0.45 -0.88 -1.12
-0.86 0.13 -1.53 0.51 2.66 1.28 -0.14 -1.19]]

代表線

上述那個圖，所有的點呈現足漸往上的趨勢。那麼，下面到底那一條線才最具資格代表所有的點呢? 藍色? 綠色 ? 紅色?

RE值

什麼是最具資格的代表線呢? 這句話太模糊了，所以需要定義一下。

定義 : 所有的點，垂直投射到某一線條時，點與線的垂直距離加總值，我們稱為RE值。若RE值為最小，則稱此線為最具代表的線，又稱為PCA線。如下圖紅色線的加總，為RE值。

上述圖形，可由如下代碼產生。

#PCA : 主成份分析 Principal Component Analysis
import numpy as np
import pylab as plt
from sklearn.decomposition import PCA

rng =np.random.RandomState(1)
W=rng.rand(2,2)#產生2列2行的亂數
#print(W)
#產生標準分佈亂數,
#scale : 愈大愈矮胖，愈小愈高瘦
X_normal=rng.normal(scale=5, size=(2,20))
X_orig=W@X_normal
X_mean=X_orig.mean(axis=1).reshape(2,1)
x=X_orig-X_mean#將整個圖型置中
#x=X_orig
print(x)
plt.figure(figsize=(6,6))
plt.xlim(-8,8)
plt.ylim(-8,8)
plt.plot([-8, 8], [0,0], linewidth=0.5, c='g')
plt.plot([0, 0], [-8,8], linewidth=0.5, c='g')
plt.scatter(x[0], x[1])

#底下 n_components = 1 表示將n維的資料降成一維(只有一項特徵)
pca_1d = PCA(n_components=1, random_state=9527)
#pca_1d.fit(x.T) #開始訓練, 取得 pca 線的向量值
#pca=pca_1d.transform(x.T).T #利用上求的向量值求每點投影到pca線時的長度
#底下是上述上個指令的簡寫
pca = pca_1d.fit_transform(x.T).T
print(pca)
#訓練完之後，才會產生pca向量值
print("components : " ,pca_1d.components_[0])
vx=pca_1d.components_[0][0]#x軸向量
vy=pca_1d.components_[0][1]#y軸向量
print((vx**2+vy**2)**0.5)#結果是1，驗証向量的斜邊是否為 1
#pca線的公式為 y=vy/vx * x
plt.plot([-8, 8], [vy/vx*-8, vy/vx*8], c='g')

#繪製垂直線
points=x.T
for i in range(len(points)):
    plt.plot([points[i][0],vx*pca[0][i]], [points[i][1], vy*pca[0][i]],c='r', linewidth=1)
plt.savefig("pca.jpg")
plt.show()

上圖中，假設綠色線的向量值為$(\begin{bmatrix}0.96913439\\0.24653303\end{bmatrix})$，則(2.89，1.52)這個點投影到紅色向量線後，座落到紅色線的長度，可以使用投影矩陣計算出來$(\begin{bmatrix}2.89, 1.52\end{bmatrix} * \begin{bmatrix}0.96913439\\0.24653303\end{bmatrix}=3.1770374623550697)$

這是一個非常神奇的事，試想想，不需要計算sin, cos耶

v = np.array([0.9691344, 0.246533])[np.newaxis, :]
L=v@X
print(L)
結果 : 
[[ 3.18 0.53 5.99 -6.53 3.81 -4.39 0.37 2.31 1.04 -4.93 -2.54 -3.28
0.46 -1.51 -3.9 -0.11 7.26 4.81 -2.2 -0.35]]

那麼，上述$(\vec{v} = \begin{bmatrix}0.96913439\\0.24653303\end{bmatrix})$是怎麼出來的呢。在Python中只需使用 PCA方法即可算出3.177這個長度，根本不需要去管他是怎麼算出來的，如下代碼所示。

pca.fit_transform(X.T) 即是以 X.T訓練PCA模型，同時得到降維後的數據，此數據就是每個點正交投射到PCA線上的長度。以上述的圖型而言，就是 3.177這個數字。

pca向量線，可以由pca_1d.components_[0]取得

pca_1d = PCA(1, random_state=9527)
pca = pca_1d.fit_transform(X.T).T
print(pca)
#繪製pca線
vx=pca_1d.components_[0][0]
vy=pca_1d.components_[0][1]
plt.plot([vx*-8, vx*8],[-vy*8, vy*8] )

結果同上
[[ 3.18 0.53 5.99 -6.53 3.81 -4.39 0.37 2.31 1.04 -4.93 -2.54 -3.28
0.46 -1.51 -3.9 -0.11 7.26 4.81 -2.2 -0.35]]

PCA向量

上面說過 pca演算法是怎麼算的不用管，直接由PCA的演算法 pca_1d.components_[0] 即可得知。但是$(\begin{bmatrix}0.96913439\\0.24653303\end{bmatrix})$這個值還是需要跟大家交代一下。此值必需使用如下三角函數及反三角函數來計算。

上述的橘色線為pca線，由(2.89,1.52)這個點投射後的長度為3.177(經由上述PCA方法算出來的)。再由X2線與pca正交後的座標為 (nx, ny)。由上述的資料可以算出 x1, 及x2。

o3角度=(o1+o2)=np.degrees(np.arcsin(1.52/x1))
o1角度=np.degrees(np.arcsin(x2/x1))

所以就可以算出o2角度=o3-o1。

即然知道o2的角度後, 那麼pca線長度為1時，(x, y)的座標值為多少呢?
x=cos(o2)
y=sin(o2)
所以 $(\vec{pca} = \begin{bmatrix}cos(o2)\\sin(o2)\end{bmatrix})$即為pca線的向量值

降維

將20個點的(x,y)二維座標(共40個數字)，投射到一條線的一維座標(變成20個數字)，稱為降維。

那麼如果是20個點的三維立體空間座標，則共有60個數字，然後投射到二維的平面時，就會變成只有40個數字。

那麼三維的立體空間可以投射成一條線嗎? 這叫廢話，當然可以。此時，60個數字就會變成20個數字了。

再進階一點，如果是20個點的100維度呢，那麼就有2000個數字，投射到一維時，還是20個數字。

Scikit-learn digits資料視覺化

底下是scikit的digits資料，共有1797筆手寫數字圖形資料，每筆資料為二維的 (8,8) = 64個像素。使用PCA降成二維後，可發現不同的數字有重疊的部份，就可以知道使用PCA去預測手寫數字會相當的不準確。

from sklearn import datasets
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt
digits = datasets.load_digits()

colors = ['black', 'blue', 'purple', 'yellow', 'white', 'red', 'lime', 'cyan', 'orange', 'gray']
pca = PCA(n_components=2)
# Fit and transform the data to the model
reduced_data_pca = pca.fit_transform(digits.data)
for i in range(len(colors)):
    x = reduced_data_pca[:, 0][digits.target == i]
    y = reduced_data_pca[:, 1][digits.target == i]
    plt.scatter(x, y, c=colors[i])
plt.legend(digits.target_names, bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.xlabel('First Principal Component')
plt.ylabel('Second Principal Component')
plt.title("PCA Scatter Plot")
plt.show()

3d圖形

from sklearn import datasets
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt
digits = datasets.load_digits()

colors = ['black', 'blue', 'purple', 'yellow', 'white', 'red', 'lime', 'cyan', 'orange', 'gray']
pca = PCA(n_components=3)
# Fit and transform the data to the model
reduced_data_pca = pca.fit_transform(digits.data)
fig=plt.figure(figsize=(10,10))
ax=fig.add_subplot(1,1,1, projection="3d")
for i in range(len(colors)):
    x = reduced_data_pca[:, 0][digits.target == i]
    y = reduced_data_pca[:, 1][digits.target == i]
    z = reduced_data_pca[:, 2][digits.target == i]
    ax.scatter(x, y, z, c=colors[i])
plt.legend(digits.target_names, bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.xlabel('First Principal Component')
plt.ylabel('Second Principal Component')
plt.title("PCA Scatter Plot")
plt.savefig("pca_3d.jpg")
plt.show()

參考 : https://leemeng.tw/essence-of-principal-component-analysis.html

主成份分析PCA