（十）降维模型

一、SVD奇异值分解（主要用于图像处理）

        SVD是一种重要的矩阵分解方法。SVD 是线代中对于实数矩阵和复数矩阵的分解，将特征分解从 半正定矩阵 推广到任意 m×n矩阵

SVD原理

代码实现：

#!/usr/bin/python
#  -*- coding:utf-8 -*-
from PIL import Image
import os
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
 
if __name__ == '__main__':
    mpl.rcParams['font.sans-serif'] = [u'simHei']
    mpl.rcParams['axes.unicode_minus'] = False
    A = Image.open('girl.jpg')
    a = np.array(A)  #转换成矩阵


#由于是彩色图像，所以3通道。a的最内层数组为三个数，分别表示RGB，用来表示一个像素
u_r, sigma_r, v_r = np.linalg.svd(a[:, :, 0])
u_g, sigma_g, v_g = np.linalg.svd(a[:, :, 1])
u_b, sigma_b, v_b = np.linalg.svd(a[:, :, 2])

plt.figure(facecolor = 'w', figsize = (10, 10))
# 奇异值个数依次取：1,2,...,12。来看看一下效果
K = 12
for k in range(1, K + 1):
        print k
        R = restore1(u_r, sigma_r, v_r, k)
        G = restore1(u_g, sigma_g, v_g, k)
        B = restore1(u_b, sigma_b, v_b, k)
        I = np.stack((R, G, B), axis = 2)
        # 将图片重构后的显示出来
        plt.subplot(3, 4, k)
        plt.imshow(I)
        plt.axis('off')
        plt.title(u'奇异值个数：%d' %  k)
 
plt.suptitle(u'SVD与图像分解', fontsize = 20)
plt.tight_layout(0.1, rect = (0, 0, 1, 0.92))
plt.show()

def restore1(u, sigma, v, k):
    m = len(u)
    n = len(v)
    a = np.zeros((m, n))
    # 重构图像
    a = np.dot(u[:, :k], np.diag(sigma[:k])).dot(v[:k, :])
    # 上述语句等价于：
    # for i in range(k):
    #     ui = u[:, i].reshape(m, 1)
    #     vi = v[i].reshape(1, n)
    #     a += sigma[i] * np.dot(ui, vi)
    a[a < 0] = 0
    a[a > 255] = 255
    return np.rint(a).astype('uint8')

二、主成分分析

主成分分析（Principal Components Analysis，PCA）是一种数据降维技术，通过正交变换将一组相关性高的变量转换为较少的彼此独立、互不相关的变量，从而减少数据的维数。

1、基本原理与思想

主成分分析是最基础数据降维方法，它只需要特征值分解，就可以对数据进行压缩、去噪，应用十分广泛。

主成分分析的目的是减少数据集变量数量，同时要保留尽可能多的特征信息；方法是通过正交变换将原始变量组转换为数量较少的彼此独立的特征变量，从而减少数据集的维数。

主成分分析方法的思想是，将高维特征（n维）映射到低维空间（k维）上，新的低维特征是在原有的高维特征基础上通过线性组合而重构的，并具有相互正交的特性，即为主成分。

通过正交变换构造彼此正交的新的特征向量，这些特征向量组成了新的特征空间。将特征向量按特征值排序后，样本数据集中所包含的全部方差，大部分就包含在前几个特征向量中，其后的特征向量所含的方差很小。因此，可以只保留前 k个特征向量，而忽略其它的特征向量，实现对数据特征的降维处理。

主成分分析方法得到的主成分变量具有几个特点：（1）每个主成分变量都是原始变量的线性组合；（2）主成分的数目大大少于原始变量的数目；（3）主成分保留了原始变量的绝大多数信息；（4）各主成分变量之间彼此相互独立。

2、算法基本流程

主成分分析的基本步骤是：对原始数据归一化处理后求协方差矩阵，再对协方差矩阵求特征向量和特征值；对特征向量按特征值大小排序后，依次选取特征向量，直到选择的特征向量的方差占比满足要求为止。

算法的基本流程如下：

（1）归一化处理，数据减去平均值；
（2）通过特征值分解，计算协方差矩阵；
（3）计算协方差矩阵的特征值和特征向量；
（4）将特征值从大到小排序；
（5）依次选取特征值最大的k个特征向量作为主成分，直到其累计方差贡献率达到要求；
（6）将原始数据映射到选取的主成分空间，得到降维后的数据。

3、代码实现

import numpy as np
from sklearn import datasets
 # 加载鸢尾花数据集中的特征作为原始数据集
    def loadIris(self):
        data = datasets.load_iris()["data"]
        return data

    # 标准化数据
    def Standard(self,data):
        # axis=0按列取均值
        mean_vector=np.mean(data,axis=0)
        return mean_vector,data - mean_vector
    def getCovMatrix(self,newData):
        # rowvar=0表示数据的每一列代表一个feature
        return np.cov(newData,rowvar=0)

    # 计算协方差矩阵的特征值和特征向量
    def getFValueAndFVector(self,covMatrix):
        fValue,fVector = np.linalg.eig(covMatrix)
        return fValue,fVector

    # 得到特征向量矩阵
    def getVectorMatrix(self,fValue,fVector,k):
        fValueSort = np.argsort(fValue)
        fValueTopN = fValueSort[:-(k + 1):-1]
        return fVector[:,fValueTopN]

    # 得到降维后的数据
    def getResult(self,data,vectorMatrix):
        return np.dot(data,vectorMatrix)
	def getFValueAndFVector(self,covMatrix):
        fValue,fVector = np.linalg.eig(covMatrix)
        return fValue,fVector
 	# 得到特征向量矩阵
   	def getVectorMatrix(self,fValue,fVector,k):
        fValueSort = np.argsort(fValue)
        fValueTopN = fValueSort[:-(k + 1):-1]
        return fVector[:,fValueTopN]
	def getResult(self,data,vectorMatrix):
        return np.dot(data,vectorMatrix)

三、因子分析

1.因子分析

（1）主要思路：降维    简化数据结构

（2）目的：将（具有错综复杂关系的）变量综合为（数量较少的）因子以再现原始变量与因子的关系，通过不同的因子，对变量进行分类消除 相关性，在信息损失最小的情况下，降维

（3）步骤

选取因子分析的变量（选相关性较大的，利于降维）――标准化处理；根据样本、估计随机向量的协方差矩阵或相关矩阵；选择一种方法――估计因子载荷阵，计算关键统计特征；进行因子旋转，使因子含义清晰化，并命名，利用因子解释变量的构成；计算每个因子在各样本上的得分，得出新的因子得分变量――进一步分析。

（4）如何分析

检验变量间偏相关度KMO值>0.6,才适合做因子分析；

调整因子个数，显示共同特征后即可命名。

代码实现：

二、实验要求

采用因子分析方法，根据48位应聘者的15项指标得分，选出6名最优秀的应聘者。

三、代码

import pandas as pd
import numpy as np
import math as math
import numpy as np
from numpy import *
from scipy.stats import bartlett
from factor_analyzer import *
import numpy.linalg as nlg
from sklearn.cluster import KMeans
from matplotlib import cm
import matplotlib.pyplot as plt
def main():
    df=pd.read_csv("./data/applicant.csv")
    # print(df)
    df2=df.copy()
    print("\n原始数据:\n",df2)
    del df2['ID']
    # print(df2)

    # 皮尔森相关系数
    df2_corr=df2.corr()
    print("\n相关系数:\n",df2_corr)

    #热力图
    cmap = cm.Blues
    # cmap = cm.hot_r
    fig=plt.figure()
    ax=fig.add_subplot(111)
    map = ax.imshow(df2_corr, interpolation='nearest', cmap=cmap, vmin=0, vmax=1)
    plt.title('correlation coefficient--headmap')
    ax.set_yticks(range(len(df2_corr.columns)))
    ax.set_yticklabels(df2_corr.columns)
    ax.set_xticks(range(len(df2_corr)))
    ax.set_xticklabels(df2_corr.columns)
    plt.colorbar(map)
    plt.show()

    # KMO测度
    def kmo(dataset_corr):
        corr_inv = np.linalg.inv(dataset_corr)
        nrow_inv_corr, ncol_inv_corr = dataset_corr.shape
        A = np.ones((nrow_inv_corr, ncol_inv_corr))
        for i in range(0, nrow_inv_corr, 1):
            for j in range(i, ncol_inv_corr, 1):
                A[i, j] = -(corr_inv[i, j]) / (math.sqrt(corr_inv[i, i] * corr_inv[j, j]))
                A[j, i] = A[i, j]
        dataset_corr = np.asarray(dataset_corr)
        kmo_num = np.sum(np.square(dataset_corr)) - np.sum(np.square(np.diagonal(A)))
        kmo_denom = kmo_num + np.sum(np.square(A)) - np.sum(np.square(np.diagonal(A)))
        kmo_value = kmo_num / kmo_denom
        return kmo_value

    print("\nKMO测度:", kmo(df2_corr))

    # 巴特利特球形检验
    df2_corr1 = df2_corr.values
    print("\n巴特利特球形检验:", bartlett(df2_corr1[0], df2_corr1[1], df2_corr1[2], df2_corr1[3], df2_corr1[4],
                                  df2_corr1[5], df2_corr1[6], df2_corr1[7], df2_corr1[8], df2_corr1[9],
                                  df2_corr1[10], df2_corr1[11], df2_corr1[12], df2_corr1[13], df2_corr1[14]))

    # 求特征值和特征向量
    eig_value, eigvector = nlg.eig(df2_corr)  # 求矩阵R的全部特征值，构成向量
    eig = pd.DataFrame()
    eig['names'] = df2_corr.columns
    eig['eig_value'] = eig_value
    eig.sort_values('eig_value', ascending=False, inplace=True)
    print("\n特征值\n：",eig)
    eig1=pd.DataFrame(eigvector)
    eig1.columns = df2_corr.columns
    eig1.index = df2_corr.columns
    print("\n特征向量\n",eig1)

    # 求公因子个数m,使用前m个特征值的比重大于85%的标准，选出了公共因子是五个
    for m in range(1, 15):
        if eig['eig_value'][:m].sum() / eig['eig_value'].sum() >= 0.85:
            print("\n公因子个数:", m)
            break

    # 因子载荷阵
    A = np.mat(np.zeros((15, 5)))
    i = 0
    j = 0
    while i < 5:
        j = 0
        while j < 15:
            A[j:, i] = sqrt(eig_value[i]) * eigvector[j, i]
            j = j + 1
        i = i + 1
    a = pd.DataFrame(A)
    a.columns = ['factor1', 'factor2', 'factor3', 'factor4', 'factor5']
    a.index = df2_corr.columns
    print("\n因子载荷阵\n", a)
    fa = FactorAnalyzer(n_factors=5)
    fa.loadings_ = a
    # print(fa.loadings_)
    print("\n特殊因子方差:\n", fa.get_communalities())  # 特殊因子方差，因子的方差贡献度 ，反映公共因子对变量的贡献
    var = fa.get_factor_variance()  # 给出贡献率
    print("\n解释的总方差（即贡献率）:\n", var)

    # 因子旋转
    rotator = Rotator()
    b = pd.DataFrame(rotator.fit_transform(fa.loadings_))
    b.columns = ['factor1', 'factor2', 'factor3', 'factor4', 'factor5']
    b.index = df2_corr.columns
    print("\n因子旋转:\n", b)

    # 因子得分
    X1 = np.mat(df2_corr)
    X1 = nlg.inv(X1)
    b = np.mat(b)
    factor_score = np.dot(X1, b)
    factor_score = pd.DataFrame(factor_score)
    factor_score.columns = ['factor1', 'factor2', 'factor3', 'factor4', 'factor5']
    factor_score.index = df2_corr.columns
    print("\n因子得分：\n", factor_score)
    fa_t_score = np.dot(np.mat(df2), np.mat(factor_score))
    print("\n应试者的五个因子得分：\n",pd.DataFrame(fa_t_score))

    # 综合得分
    wei = [[0.50092], [0.137087], [0.097055], [0.079860], [0.049277]]
    fa_t_score = np.dot(fa_t_score, wei) / 0.864198
    fa_t_score = pd.DataFrame(fa_t_score)
    fa_t_score.columns = ['综合得分']
    fa_t_score.insert(0, 'ID', range(1, 49))
    print("\n综合得分：\n", fa_t_score)
    print("\n综合得分：\n", fa_t_score.sort_values(by='综合得分', ascending=False).head(6))

    plt.figure()
    ax1=plt.subplot(111)
    X=fa_t_score['ID']
    Y=fa_t_score['综合得分']
    plt.bar(X,Y,color="#87CEFA")
    # plt.bar(X, Y, color="red")
    plt.title('result00')
    ax1.set_xticks(range(len(fa_t_score)))
    ax1.set_xticklabels(fa_t_score.index)
    plt.show()

    fa_t_score1=pd.DataFrame()
    fa_t_score1=fa_t_score.sort_values(by='综合得分',ascending=False).head()
    ax2 = plt.subplot(111)
    X1 = fa_t_score1['ID']
    Y1 = fa_t_score1['综合得分']
    plt.bar(X1, Y1, color="#87CEFA")
    # plt.bar(X1, Y1, color='red')
    plt.title('result01')
    plt.show()

if __name__ == '__main__':
    main()

 

四、实验步骤

（1）引入数据，数据标准化

因为数据是面试中的得分，量纲相同，并且数据的分布无异常值，所以数据可以不进行标准化。

（2）建立相关系数矩阵

计算皮尔森相关系数，从热图中可以明显看出变量间存在的相关性。

 

进行相关系数矩阵检验——KMO测度和巴特利特球体检验：

KMO值：0.9以上非常好；0.8以上好；0.7一般；0.6差；0.5很差；0.5以下不能接受；巴特利球形检验的值范围在0-1，越接近1，使用因子分析效果越好。

通过观察上面的计算结果，可以知道，KMO值为0.783775605643526，在较好的范围内，并且巴特利球形检验的值接近1，所有可以使用因子分析。

（3）求解特征值及相应特征向量 

 

求公因子个数m,使用前m个特征值的比重大于85%的标准，选出了公共因子是五个。

（4）因子载荷阵

  

由上可以看出，选择5个公共因子，从方差贡献率可以看出，其中第一个公因子解释了总体方差的50.092%，四个公共因子的方差贡献率为86.42%，可以较好的解释总体方差。

（5）因子旋转

 

（6）因子得分

 

（7）根据应聘者的五个因子得分，按照贡献率进行加权，得到最终各应试者的综合得分，然后选出前六个得分最高的应聘者。