本文档是对课程内容的简单整理。

    • 前缀Z_ 为正课视频的简称，如Z1为正课视频第1课
    • 前缀G_ 为更新视频的简称，如G1为更新视频第1课

    主成分分析
    Matlab 见D:\00000MCM\清风\0 课件和代码\正课配套的课件和代码\第14讲.主成分分析\代码和例题数据
    SPSS 因子分析（不进行旋转）数学建模——主成分分析(SPSS） - 知乎
    Python 见D:\00000py\1\mcm 2019d_q3_0030_PCAKmeans.py

    主要使用Matlab or SPSS

    文章 D:\00000MCM\0 研究生数学建模竞赛历年真题和优秀论文集锦\研究生数学建模-优秀论文\2019优秀数模论文\D19102470244.pdf（用于降维）
    0244 0030
    2018C 1 PCA 用于降维

    Python

    1. from sklearn.decomposition import PCA
    2. from sklearn.cluster import KMeans
    3. import matplotlib.pyplot as plt
    4. from mpl_toolkits.mplot3d import Axes3D
    5. import numpy as np
    6. import pandas as pd
    7. import copy
    8. import datetime
    10. paths = './datasets/2019d/'
    12. data1 = pd.read_excel(paths + 'feature1_forq3.xlsx', header=None)    #文件 1 特征参数数据
    13. data1 = pd.DataFrame(data1)
    14. data1 = np.array(data1)
    15. data0 = pd.read_excel(paths + 'data1_small_outforq2.xlsx', header=None)    #文件 1汽车运行数据
    16. data0 = pd.DataFrame(data0)
    17. data0 = np.array(data0)
    18. sport = pd.read_excel(paths + 'sport1_forq3.xlsx', header=None)  # 文件 1 运动学片段
    19. sport = pd.DataFrame(sport)
    20. sport = np.array(sport)
    21. sport = sport - np.array([[1, 1, 0]] * sport.shape[0])
    23. # 合并数据 转换数据
    24. data = data1.tolist()
    25. data = np.array(data)  # (200, 15)   200行数据 15特征
    27. #PCA 3 个主成分
    28. pca_d = PCA(n_components=3)    #选择 3 个主成分
    29. data_new = pca_d.fit_transform(data)  # (200, 3)  200行数据 3主成分
    30. clust = KMeans(n_clusters=3)    #构造聚类,分 3 类
    31. clust.fit(data_new)    #聚类
    32. label_p = clust.labels_    #获取聚类标签  (200, )  200行数据 1标签     范围 0 1 2
    34. la1 = data_new[label_p == 0]   # (170, 3)
    35. la2 = data_new[label_p == 1]   # (24, 3)
    36. la3 = data_new[label_p == 2]   # (6, 3)
    37. x1=[]
    38. y1=[]
    39. z1=[]
    40. x2=[]
    41. y2=[]
    42. z2=[]
    43. x3=[]
    44. y3=[]
    45. z3=[]
    46. for i in range(la1.shape[0]):  # 三位空间分布
    47.     x1.append(la1[i][0])
    48.     y1.append(la1[i][1])
    49.     z1.append(la1[i][2])
    50. for i in range(la2.shape[0]):
    51.     x2.append(la2[i][0])
    52.     y2.append(la2[i][1])
    53.     z2.append(la2[i][2])
    54. for i in range(la3.shape[0]):
    55.     x3.append(la3[i][0])
    56.     y3.append(la3[i][1])
    57.     z3.append(la3[i][2])
    58. fig = plt.figure()
    59. ax = fig.add_subplot(111, projection='3d')
    60. ax.scatter(x1,y1,z1,c='b')
    61. ax.scatter(x2,y2,z2,c='r')
    62. ax.scatter(x3,y3,z3,c='k')
    63. plt.show()

    Matlab

    1. clear;clc
    2. load data1.mat   % 主成分聚类
    3. %  load data2.mat   % 主成分回归
    5. % 注意，这里可以对数据先进行描述性统计
    6. % 描述性统计的内容见第5讲.相关系数
    7. [n,p] = size(x);  % n是样本个数，p是指标个数     n行样本个数 p列特征
    9. %% 第一步：对数据x标准化为X
    10. X=zscore(x);   % matlab内置的标准化函数（x-mean(x)）/std(x)
    12. %% 第二步：计算样本协方差矩阵
    13. R = cov(X);
    15. %% 注意：以上两步可合并为下面一步：直接计算样本相关系数矩阵
    16. R = corrcoef(x);
    17. disp('样本相关系数矩阵为：')
    18. disp(R)
    20. %% 第三步：计算R的特征值和特征向量
    21. % 注意：R是半正定矩阵，所以其特征值不为负数
    22. % R同时是对称矩阵，Matlab计算对称矩阵时，会将特征值按照从小到大排列哦
    23. % eig函数的详解见第一讲层次分析法的视频
    24. [V,D] = eig(R);  % V 特征向量矩阵  D 特征值构成的对角矩阵
    27. %% 第四步：计算主成分贡献率和累计贡献率
    28. lambda = diag(D);  % diag函数用于得到一个矩阵的主对角线元素值(返回的是列向量)
    29. lambda = lambda(end:-1:1);  % 因为lambda向量是从小大到排序的，我们将其调个头
    30. contribution_rate = lambda / sum(lambda);  % 计算贡献率
    31. cum_contribution_rate = cumsum(lambda)/ sum(lambda);   % 计算累计贡献率  cumsum是求累加值的函数
    32. disp('特征值为：')
    33. disp(lambda')  % 转置为行向量，方便展示
    34. disp('贡献率为：')
    35. disp(contribution_rate')
    36. disp('累计贡献率为：')
    37. disp(cum_contribution_rate')
    38. disp('与特征值对应的特征向量矩阵为：')
    39. % 注意：这里的特征向量要和特征值一一对应，之前特征值相当于颠倒过来了，因此特征向量的各列需要颠倒过来
    40. %  rot90函数可以使一个矩阵逆时针旋转90度，然后再转置，就可以实现将矩阵的列颠倒的效果
    41. V=rot90(V)';
    42. disp(V)
    45. %% 计算我们所需要的主成分的值
    46. m =input('请输入需要保存的主成分的个数:  ');
    47. F = zeros(n,m);  %初始化保存主成分的矩阵（每一列是一个主成分）
    48. for i = 1:m
    49.     ai = V(:,i)';   % 将第i个特征向量取出，并转置为行向量
    50.     Ai = repmat(ai,n,1);   % 将这个行向量重复n次，构成一个n*p的矩阵
    51.     F(:, i) = sum(Ai .* X, 2);  % 注意，对标准化的数据求了权重后要计算每一行的和
    52. end
    54. %% (1)主成分聚类 ： 将主成分指标所在的F矩阵复制到Excel表格，然后再用Spss进行聚类
    55. % 在Excel第一行输入指标名称（F1,F2, ..., Fm）
    56. % 双击Matlab工作区的F,进入变量编辑中，然后复制里面的数据到Excel表格
    57. % 导出数据之后，我们后续的分析就可以在Spss中进行。
    59. %%（2）主成分回归：将x使用主成分得到主成分指标，并将y标准化，接着导出到Excel，然后再使用Stata回归
    60. % Y = zscore(y);  % 一定要将y进行标准化哦~
    61. % 在Excel第一行输入指标名称（Y,F1, F2, ..., Fm）
    62. % 分别双击Matlab工作区的Y和F,进入变量编辑中，然后复制里面的数据到Excel表格
    63. % 导出数据之后，我们后续的分析就可以在Stata中进行。

    image.png

    降维
    主成分分析用于聚类 - 可视化
    主成分回归