定义
输入: x = ( x 1 , x 2 , ⋯ , x m ) T , x : m 维随机变量，且满足均值向量 μ = E ( x ) = ( μ 1 , μ 2 , ⋯ , μ m ) T , 协方差矩阵 ∑ = c o v ( x , x ) = E [ ( x − μ ) ( x − μ ) T ] x = (x_1,x_2,\cdots,x_m)^T,x:m维随机变量，且满足均值向量\mu = E(x) = (\mu_1,\mu_2,\cdots,\mu_m)^T,协方差矩阵\sum = cov(x,x) = E\big[ (x - \mu)(x - \mu)^T \big] x=(x1,x2,⋯,xm)T,x:m维随机变量，且满足均值向量μ=E(x)=(μ1,μ2,⋯,μm)T,协方差矩阵∑=cov(x,x)=E[(x−μ)(x−μ)T]
输出: y = ( y 1 , y 2 , ⋯ , y m ) T , y i = α i T x = α 1 i x 1 + α 2 i x 2 + ⋯ + α m i x m , α i T = ( α 1 i , α 2 i , ⋯ , α m i ) , i = 1 , 2 , ⋯ , m 。 y = (y_1,y_2,\cdots,y_m)^T,y_i=\alpha_i^Tx=\alpha_{1i}x_1+\alpha_{2i}x_2+\cdots+\alpha_{mi}x_m,\alpha_i^T = (\alpha_{1i},\alpha_{2i},\cdots,\alpha_{mi}),i=1,2,\cdots,m。 y=(y1,y2,⋯,ym)T,yi=αiTx=α1ix1+α2ix2+⋯+αmixm,αiT=(α1i,α2i,⋯,αmi),i=1,2,⋯,m。
(1)系数向量 α i T \alpha_i^T αiT是单位向量,即 α i T α i = 1 , i = 1 , 2 , ⋯ , m ; \alpha_i^T\alpha_i = 1,i=1,2,\cdots,m; αiTαi=1,i=1,2,⋯,m;
(2)变量 y i y_i yi与 y j y_j yj互不相关,即 c o v ( y i , y j ) = 0 ( i ≠ j ) ; cov(y_i,y_j) = 0(i\ne j); cov(yi,yj)=0(i=j);
(3)变量 y 1 y_1 y1是 x x x的所有线性变换中方差最大的; y 2 y_2 y2是与 y 1 y_1 y1不相关的 x x x的所有线性变换中方差最大的;一般地, y i y_i yi是与 y 1 , y 2 , ⋯ , y i − 1 ( i = 1 , 2 , ⋯ , m ) y_1,y2,\cdots,y_{i-1}(i=1,2,\cdots,m) y1,y2,⋯,yi−1(i=1,2,⋯,m)都不相关的 x x x的所有线性变换中方差最大的;这时分别称 y 1 , y 2 , ⋯ , y m y1,y_2,\cdots,y_m y1,y2,⋯,ym为 x x x的 第一主成分、第二主成分、 ⋯ 、第 m 主成分 第一主成分、第二主成分、\cdots、第m主成分 第一主成分、第二主成分、⋯、第m主成分。
输入空间
T= { ( x 1 , x 2 , , … , x N ) } \left\{(x_1,x_2,,\dots,x_N)\right\} {(x1,x2,,…,xN)}
import numpy as np
import pandas as pd
import time def load_data(file):'''数据集 下载地址：https://download.csdn.net/download/nanxiaotao/89743722INPUT:file - (str) 数据文件的路径 OUTPUT:df - (dataframe) 读取的数据表格X - (array) 特征数据数组'''df = pd.read_csv(file) #读取csv文件df.drop('Sports', axis=1, inplace=True) #去掉类别数据X = np.asarray(df.values).T #将数据转换成数组return df, X
df, X = load_data('cars.csv') #加载数据
np.shape(X)
def Normalize(X):'''规范化INPUT:X - (array) 特征数据数组OUTPUT:X - (array) 规范化处理后的特征数据数组'''m, n = X.shapefor i in range(m):E_xi = np.mean(X[i]) #第i列特征的期望Var_xi = np.var(X[i], ddof=1) #第i列特征的方差#print(i,np.isnan(Var_xi))for j in range(n):X[i][j] = (X[i][j] - E_xi) / np.sqrt(Var_xi) #对第i列特征的第j条数据进行规范化处理return X
X = Normalize(X) #对样本数据进行规范化处理
np.shape(X)
特征空间（Feature Space）
X[0:16,0]
算法
A = U ∗ ∑ V T , A ∈ R m x n , U : m 阶正交矩阵 , V : n 阶正交矩阵 , ∑ : 由降序排列的非负的对角线元素组成的 m ∗ n 矩形对角矩阵 A = U* \sum V^T,A \in R^{mxn},U:m阶正交矩阵,V:n阶正交矩阵,\sum:由降序排列的非负的对角线元素组成的m*n矩形对角矩阵 A=U∗∑VT,A∈Rmxn,U:m阶正交矩阵,V:n阶正交矩阵,∑:由降序排列的非负的对角线元素组成的m∗n矩形对角矩阵
U U T = I , V V T = I , I : 单位矩阵 UU^T = I,VV^T = I,I:单位矩阵 UUT=I,VVT=I,I:单位矩阵
∑ = d i a g ( σ 1 , σ 2 , ⋯ , σ p ) , σ 1 ≥ σ 2 ≥ ⋯ ≥ σ p ≥ 0 , p = m i n ( m , n ) \sum = diag(\sigma_1,\sigma_2,\cdots,\sigma_p),\sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_p \geq 0,p = min(m,n) ∑=diag(σ1,σ2,⋯,σp),σ1≥σ2≥⋯≥σp≥0,p=min(m,n)
def cal_V(X):'''奇异值分解：计算V矩阵和特征值INPUT:X - (array) 特征数据数组OUTPUT:eigvalues - (list) 特征值列表，其中特征值按从大到小排列V - (array) V矩阵'''newX = X.T / np.sqrt(X.shape[1]-1) #构造新矩阵X'Sx = np.matmul(newX.T, newX) #计算X的协方差矩阵Sx = X'.T * X'V_T = [] #用于保存V的转置#print(np.isinf(Sx))#print(np.isnan(Sx))w, v = np.linalg.eig(Sx) #计算Sx的特征值和对应的特征向量，即为X’的奇异值和奇异向量tmp = {} #定义一个字典用于保存特征值和特征向量，字典的键为特征值，值为对应的特征向量for i in range(len(w)):tmp[w[i]] = v[i]eigvalues = sorted(tmp, reverse=True) #将特征值逆序排列后保存到eigvalues列表中for i in eigvalues:d = 0for j in range(len(tmp[i])):d += tmp[i][j] ** 2V_T.append(tmp[i] / np.sqrt(d)) #计算特征值i的单位特征向量，即为V矩阵的列向量，将其保存到V_T中V = np.array(V_T).T #对V_T进行转置得到V矩阵return eigvalues, V
y k = α k T x = α 1 k x 1 + α 2 k x 2 + ⋯ + α m k x m , α k T = ( α 1 k , α 2 k , ⋯ , α m k ) , k = 1 , 2 , ⋯ , m y_k=\alpha_k^Tx=\alpha_{1k}x_1+\alpha_{2k}x_2+\cdots+\alpha_{mk}x_m,\alpha_k^T = (\alpha_{1k},\alpha_{2k},\cdots,\alpha_{mk}),k=1,2,\cdots,m yk=αkTx=α1kx1+α2kx2+⋯+αmkxm,αkT=(α1k,α2k,⋯,αmk),k=1,2,⋯,m
v a r ( y k ) = α k T ∑ α k = λ k , k = 1 , 2 , ⋯ , m var(y_k) = \alpha_k^T\sum \alpha_k = \lambda_k,k = 1,2,\cdots,m var(yk)=αkT∑αk=λk,k=1,2,⋯,m
m a x α 1 \mathop{max}\limits_{\alpha_1} α1max α 1 T ∑ α 1 \alpha_1^T\sum\alpha_1 α1T∑α1
s . t . s.t. s.t. α 1 T α 1 = 1 \alpha_1^T\alpha_1 = 1 α1Tα1=1
def do_pca(X, k):'''主成分分析INPUT:X - (array) 特征数据数组k - (int) 设定的主成分个数OUTPUT:fac_load - (array) 因子负荷量数组dimrates - (list) 可解释偏差列表Y - (array) 主成分矩阵'''eigvalues, V = cal_V(X)Vk = V[:, :k] Y = np.matmul(Vk.T, X) dimrates = [i / sum(eigvalues) for i in eigvalues[:k]]fac_load = np.zeros((k, X.shape[0]))for i in range(k): for j in range(X.shape[0]):fac_load[i][j] = np.sqrt(eigvalues[i]) * Vk[j][i] / np.sqrt(np.var(X[j]))return fac_load, dimrates, Y
k = 3 #设定主成分个数为3
fac_load, dimrates, Y = do_pca(X, k) #进行主成分分析
pca_result = pd.DataFrame(fac_load, index=['Dimension1', 'Dimension2', 'Dimension3'], columns=df.columns) #将结果保存为dataframe格式
pca_result['Explained Variance'] = dimrates #将可解释偏差保存到pca_result的'Explained Variance'列
pca_result