数学建模-主成分析

在实际问题研究中，多变量问题是经常会遇到的。变量太多，无疑会增加分析问题的难度与复杂性，而且在许多实际问题中，多个变量之间是具有一定的相关关系的。 因此，人们会很自然地想到，能否在相关分析的基础上，用较少的新变量代替原来较多的旧变量，而且使这些较少的新变量尽可能多地保留原来变量所反映的信息。

主成分分析是把原来多个变量划为少数几个综合指标的一种统计分析方法。

主成分析核心思想：降维

一般取累计贡献率达85%—95%的特征值对应的向量作为主成分

建模过程：

(8条消息) PCA：详细解释主成分分析_lanyuelvyun的博客-CSDN博客_pca

数学建模——主成分分析(SPSS） - 知乎 (zhihu.com)

在建模中主要有两种方法可以实现主成分析：

1、MATLAB

clear;clc
x = xlsread('主成分分析_学生成绩样例.xlsx',1,'E5:M14')  % 读数据

[n,p] = size(x);  % n是样本个数，p是指标个数
                  % n返回的是行数，p返回的是列数
                  
% 第一步：数据x标准化
X=zscore(x);      % zscore()：matlab内置的标准化函数（x-mean(x)）/std(x)

% 第二步：计算样本协方差矩阵
C = cov(X);

% 其实以上两步可合并为下面的一步：直接计算样本相关系数矩阵
C = corrcoef(x);

disp('样本相关系数矩阵为：')
disp(C);

% 第三步：计算C的特征值和特征向量

% 注意：C是半正定矩阵，所以其特征值不为负数
% C同时是对称矩阵，Matlab在处理对称矩阵时，会将特征值按照从小到大排列
[V,D] = eig(C);  % eig()：计算特征值和特征向量
                 % V 特征向量矩阵  D是由特征值构成的对角矩阵
                 

% 第四步：计算主成分贡献率和累计贡献率
lambda = diag(D);           % diag()：用于得到一个矩阵的主对角线元素值(返回的是列向量)

lambda = lambda(end:-1:1);  % 因为lambda向量是从小大到排序的，我们将其调个顺序

contribution_rate = lambda / sum(lambda);              % 计算贡献率

cum_contribution_rate = cumsum(lambda)/ sum(lambda);   % 计算累计贡献率  cumsum()：求累加值的函数

disp('特征值为：')
disp(lambda')        


disp('贡献率为：')
disp(contribution_rate')


disp('累计贡献率为：')
disp(cum_contribution_rate')


disp('与特征值对应的特征向量矩阵为：')   % 注意：这里的特征向量要和特征值一一对应，之前特征值相当于颠倒过来了，因此特征向量的各列需要颠倒过来
V=rot90(V)';                          % rot90()：可以使一个矩阵逆时针旋转90度，然后再转置，就可以实现将矩阵的列颠倒的效果
disp(V)

xlswrite('主成分分析_学生成绩样例.xlsx',lambda,2);   % 把特征值写入当前文件的 脚本2 里面
xlswrite('主成分分析_学生成绩样例.xlsx',V,3);        % 把特征向量写入当前文件的 脚本3 里面

2、SPSS

(8条消息) 数学建模：评价模型——主成分分析 PCA SPSS实现、python实现_Luweir23的博客-CSDN博客_主成分分析评价指标

3、Python

1、简单主成分析

import numpy as np
from sklearn.decomposition import PCA

a = np.loadtxt("Pdata11_7.txt")
b = np.r_[a[:, 1:4], a[:, -3:]]  # 构造数据矩阵
print("相关系数矩阵：", np.around(np.corrcoef(b.T), decimals=3))  # 数据标准化并计算相关系数矩阵，并保留三位小数

md = PCA(n_components=0.85).fit(b)  # 构造并训练模型(累计贡献率>85%即可)

print("特征值为：", md.explained_variance_)
print("各主成分的贡献率：", md.explained_variance_ratio_)
print("奇异值为：", md.singular_values_)
print("各主成分的系数：\n", md.components_)  # 每行是一个主成分
"""下面直接计算特征值和特征向量，和库函数进行对比"""
cf = np.cov(b.T)  # 计算协方差阵
c, d = np.linalg.eig(cf)  # 求特征值和特征向量
print("特征值为：", c)
print("特征向量为：\n", d)
print("各主成分的贡献率为：", c / np.sum(c))

2、综合评价

import numpy as np
from scipy.stats import zscore

a = np.loadtxt("Pdata11_8.txt")
print("相关系数阵为：\n", np.corrcoef(a.T))
b = np.delete(a, 0, axis=1)  # 删除第1列数据
c = zscore(b)
r = np.corrcoef(c.T)  # 数据标准化并计算相关系数阵
d, e = np.linalg.eig(r)  # 求特征值和特征向量
rate = d / d.sum()  # 计算各主成分的贡献率
print("特征值为：", d)
print("特征向量为：\n", e)
print("各主成分的贡献率为：", rate)
k = 1  # 提出主成分的个数
F = e[:, :k]
score_mat = c.dot(F)  # 计算主成分得分矩阵
score1 = score_mat.dot(rate[0:k])  # 计算各评价对象的得分
score2 = -score1  # 通过表中数据以及score1观测，需要调整得分的正负号
print("各评价对象的得分为：", score2)
index = score1.argsort() + 1  # 排序后的每个元素在原数组中的位置
print("从高到低各个城市的编号排序为：", index)