• 数学建模--数据预处理

    目录

    一、数据统计

    1、行列式的最大元素和最小元素

    2、求向量的平均值和中值 

    3、对矩阵做排序

     二、数据优化（数据残缺值和异常值的处理）

    1、数据残缺

    ①插值

    ②拟合

     ③邻近替换

    ④KNN算法填充

    2、数据异常

    ①拉依达准则

    ②替换异常值

    3、数据变换

    ①0-1标准化

    ②z-score标准化

    ③标准化的应用

    4、数据离散化

    ①等宽法

    ②等频法

    三、数据降维

    1、主成分分析

    ①主成分分析简介

    ②主成分分析计算步骤

    ③案例分析

     2、因子分析（FA）

    ①应用场景

    ②步骤

    ③模型分析

     ④代码实现

    一、数据统计

    数据统计一般包括求矩阵最大、最小元素，求矩阵平均值和中值， 矩阵元素求和、求积，矩阵元素累加和与累乘积，求标准方差、相关系数、元素排序等。

    直接举例子说明

    1、行列式的最大元素和最小元素

    1. 命令如下：
    2. A=[12,45,58;25,60,-45;56,25,178;2,0,-13];
    3. max(A,[],2)        %求每行最大元素
    4. ans =
    5.         58
    6.         60
    7.         178
    8.         2
    9. max(A)             %求每列最大元素
    10. ans =
    11. 56    60   178
    12. min(min(A))      　　 %求整个矩阵的最小元素。也可用命令：min(A(：))
    13. ans =
    14. -45

    2、求向量的平均值和中值 

    1. 命令如下:
    2. x=[10,24,6,-8,0,-12];         %偶数个元素
    3. mean(x)                   %求此向量的平均值
    4. x =
    5. 10    24     6    -8     0   -12
    6. ans =
    7. 3.3333
    8. median(x)                  %求此向量的中值
    9. ans =
    10. 3

    3、对矩阵做排序

    1. 命令如下：
    2.  
    3.     A=[0,-11,5;1,15,7;16,9,-20];
    4.  
    5. sort(A)                    %对A的每列按升序排序
    6.  
    7. ans =
    8.  
    9.         0    -11  -20
    10.  
    11.         1     9    5
    12.  
    13.        16    15    7
    14.  
    15.  
    16.  
    17. -sort(-A,2)                 %对A的每行按降序排列
    18.  
    19. ans =
    20.  
    21.         5     0   -11
    22.  
    23.        15     7     1
    24.  
    25.        16     9   -20
    26.  
    27.  
    28.  
    29. [X,I]=sort(A)               %对A按列排序，并将每个元素所在行号送矩阵I
    30.  
    31. X =
    32.  
    33.        0    -11   -20
    34.  
    35.        1     9     5
    36.  
    37.       16    15     7
    38.  
    39. I =
    40.  
    41.        1     1     3
    42.  
    43.        2     3     1
    44.  
    45.        3     2     2

     二、数据优化（数据残缺值和异常值的处理）

    插值和拟合都是数据优化的一种方法，当实验数据不够多时经常需要用到这种方法来画图。在MATLAB中都有特定的函数来完成这些功能。这两种方法的确别在于：

    当测量值是准确的，没有误差时，一般用插值；

    当测量值与真实值有误差时，一般用数据拟合。

    1、数据残缺

    ①插值

    对于一维曲线的插值，一般用到的函数yi=interp1(X,Y,xi,method) ，其中method包括nearst，linear，spline，cubic。

    对于二维曲面的插值，一般用到的函数zi=interp2(X,Y,Z,xi,yi,method)，其中method也和上面一样，常用的是cubic。

    1. %产生原始数据
    2. x=0:0.1:1;
    3. y=(x.^2-3*x+7).*exp(-4*x).*sin(2*x);
    4. subplot(2,2,1);
    5. plot(x,y);
    6. title('原始数据');
    7. %线性插值
    8. xx=0:0.01:1;
    9. y1=interp1(x,y,xx,'linear');
    10. %subplot(2,2,1)
    11. %plot(x,y,'o',xx,y1);
    12. %title('线性插值');
    13. %最邻近点插值
    14. y2=interp1(x,y,xx,'nearest');
    15. subplot(2,2,2)
    16. plot(x,y,'o',xx,y2);
    17. title('最邻近点插值');
    18. %三次插值
    19. y3=interp1(x,y,xx,'pchip');
    20. subplot(2,2,3)
    21. plot(x,y,'o',xx,y3);
    22. title('三次插值');
    23. %三次样条插值
    24. y4=interp1(x,y,xx,'spline');
    25. subplot(2,2,4)
    26. plot(x,y,'o',xx,y4);
    27. title('三次样条插值');

    经典问题：利用给定的高度补充地图

    1. %插值基点为网格节点
    2. clear all
    3. y=20:-1:0;
    4. x=0:20;
    5. z=[0.2	0.2	0.2	0.2	0.2	0.2	0.4	0.4	0.3	0.2	0.3	0.2	0.1	0.2	0.2	0.4	0.3	0.2	0.2	0.2	0.2;
    6. 0.3	0.2	0.2	0.2	0.2	0.4	0.3	0.3	0.3	0.3	0.4	0.2	0.2	0.2	0.2	0.4	0.4	0.4	0.3	0.2	0.2;
    7. 0.2	0.3	0.3	0.2	0.3	1	0.4	0.5	0.3	0.3	0.3	0.3	0.2	0.2	0.2	0.6	0.5	0.4	0.4	0.2	0.2;
    8. 0.2	0.2	0.4	0.2	1	1.1	0.9	0.4	0.3	0.3	0.5	0.3	0.2	0.2	0.2	0.7	0.3	0.6	0.6	0.3	0.4;
    9. 0.2	0.2	0.9	0.7	1	1	1	0.7	0.5	0.3	0.2	0.2	0.2	0.6	0.2	0.8	0.7	0.9	0.5	0.5	0.4;
    10. 0.2	0.3	1	1	1	1.2	1	1.1	0.8	0.3	0.2	0.2	0.2	0.5	0.3	0.6	0.6	0.8	0.7	0.6	0.5;
    11. 0.2	0.4	1	1	1.1	1.1	1.1	1.1	0.6	0.3	0.4	0.4	0.2	0.7	0.5	0.9	0.7	0.4	0.9	0.8	0.3;
    12. 0.2	0.2	0.9	1.1	1.2	1.2	1.1	1.1	0.6	0.3	0.5	0.3	0.2	0.4	0.3	0.7	1	0.7	1.2	0.8	0.4;
    13. 0.2	0.3	0.4	0.9	1.1	1	1.1	1.1	0.7	0.4	0.4	0.4	0.3	0.5	0.5	0.8	1.1	0.8	1.1	0.9	0.3;
    14. 0.3	0.3	0.5	1.2	1.2	1.1	1	1.2	0.9	0.5	0.6	0.4	0.6	0.6	0.3	0.6	1.2	0.8	1	0.8	0.5;
    15. 0.3	0.5	0.9	1.1	1.1	1	1.2	1	0.8	0.7	0.5	0.6	0.4	0.5	0.4	1	1.3	0.9	0.9	1	0.8;
    16. 0.3	0.5	0.6	1.1	1.2	1	1	1.1	0.9	0.4	0.4	0.5	0.5	0.8	0.6	0.9	1	0.5	0.8	0.8	0.9;
    17. 0.4	0.5	0.4	1	1.1	1.2	1	0.9	0.7	0.5	0.6	0.3	0.6	0.4	0.6	1	1	0.6	0.9	1	0.7;
    18. 0.3	0.5	0.8	1.1	1.1	1	0.8	0.7	0.7	0.4	0.5	0.4	0.4	0.5	0.4	1.1	1.3	0.7	1	0.7	0.6;
    19. 0.3	0.5	0.9	1.1	1	0.7	0.7	0.4	0.6	0.4	0.4	0.3	0.5	0.5	0.3	0.9	1.2	0.8	1	0.8	0.4;
    20. 0.2	0.3	0.6	0.9	0.8	0.8	0.6	0.3	0.4	0.5	0.4	0.5	0.4	0.2	0.5	0.5	1.3	0.6	1	0.9	0.3;
    21. 0.2	0.3	0.3	0.7	0.6	0.6	0.4	0.2	0.3	0.5	0.8	0.8	0.3	0.2	0.2	0.8	1.3	0.9	0.8	0.8	0.4;
    22. 0.2	0.3	0.3	0.6	0.3	0.4	0.3	0.2	0.2	0.3	0.6	0.4	0.3	0.2	0.4	0.3	0.8	0.6	0.7	0.4	0.4;
    23. 0.2	0.3	0.4	0.4	0.2	0.2	0.2	0.3	0.2	0.2	0.2	0.2	0.2	0.2	0.2	0.5	0.7	0.4	0.4	0.3	0.3;
    24. 0.2	0.2	0.3	0.2	0.2	0.3	0.2	0.2	0.2	0.2	0.2	0.1	0.2	0.4	0.3	0.6	0.5	0.3	0.3	0.3	0.2;
    25. 0.2	0.2	0.2	0.2	0.2	0.2	0.2	0.2	0.2	0.2	0.2	0.2	0.2	0.2	0.4	0.7	0.4	0.2	0.4	0.5	0.5];
    26.  
    27. %未插值直接画图
    28. figure(1)                                         %创建图形窗口1，并激活
    29. surf(x,y,z);
    30. shading flat                                      %用shading flat命令，使曲面变的光滑
    31. title('未插值地形图')
    32. xlabel('横坐标')
    33. ylabel('纵坐标')
    34. zlabel('高度')
    35.  
    36. %三次插值后画图
    37. %画地形图
    38. figure(2)
    39. xi=0:0.05:20;
    40. yi=20:-0.05:0;
    41. zi=interp2(x,y,z,xi',yi,'cubic');                 %'cubic'三次插值
    42. surfc(xi,yi,zi);                                  %底面带等高线
    43. shading flat
    44. title('插值后地形图')
    45. xlabel('横坐标')
    46. ylabel('纵坐标')
    47. zlabel('高度')
    49. %画立体等高线图
    50. figure(3)
    51. contour3(xi,yi,zi);
    52. title('立体等高线图')
    53. xlabel('横坐标')
    54. ylabel('纵坐标')
    55. zlabel('高度')
    57. %画等高线图
    58. figure(4)
    59. [c,h]=contour(xi,yi,zi);
    60. clabel(c,h);                                       %用于为2维等高线添加标签
    61. colormap cool                                     %冷色调
    62. title('平面等高线图')
    63. xlabel('横坐标')
    64. ylabel('纵坐标')

    ②拟合

    对于一维曲线的拟合，一般用到的函数p=polyfit(x,y,n)和yi=polyval(p,xi)，这个是最常用的最小二乘法的拟合方法。

    对于二维曲面的拟合，有很多方法可以实现，这里运用Spline Toolbox里面的函数功能。

    1. x = 0:0.1:1;
    2. y = [-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2];
    3. A = polyfit(x,y,2)   %A为拟合出来的函数
    4. z = polyval(A,x);    %求多项式在x处的值z
    5. plot(x,y,'k+',x,z,'r')

     ③邻近替换

    前/后一个非缺失值将其替换

    最近的非缺失值替换

    相邻的非离群值线性插值填充

    1. test_data1=fillmissing(test_data,'previous');
    2. test_data1=fillmissing(test_data,'next');
    3. test_data1=fillmissing(test_data,'nearest');
    4. test_data1=fillmissing(test_data,'linear');

    ④KNN算法填充

    1. from sklearn.metrics import nan_euclidean_distances
    2. import numpy as np
    3. from sklearn.impute import KNNImputer
    4.  
    5. X = [[1, 2, np.nan], [3, 4, 3], [np.nan, 6, 5], [8, 8, 7]]
    6.  
    7. # python的nan_euclidean_distances函数可计算含空值的距离矩阵
    8. nan_euclidean_distances(X, X)
    9.  
    10. # python用KNNImputer进行空值填充
    11. X = [[1, 2, np.nan], [3, 4, 3], [np.nan, 6, 5], [8, 8, 7]]
    12. imputer = KNNImputer(n_neighbors=2)
    13. a = imputer.fit_transform(X)
    14. print(a)

    2、数据异常

    ①拉依达准则

    拉依达准则是用来发现数据异常值

    1. x=[1, 1.1, 1.2, 1.3, 1.4, 2, 1.2, 1.3, 1.5, 0.9, 0.8, 1.1, 11];
    2. inlier = [];outlier = [];
    3. len = length(x);
    4. average1 = mean(x);  % x中所有元素的均值
    5. standard1 = std(x);  % x的标准差
    6. for i = 1:len  % 遍历x向量，判断是否为偏离点，不是偏离点则存入inline
    7.      if abs(x(i)-average1)3
    8.        inlier = [inlier x(i)];
    9.      end
    10.    end
    11.  
    12.  
    13. average2 = mean(inlier);
    14. standard2 = std(inlier);
    15. for i = 1:len  % 遍历x向量，判断是否为偏离点，不是偏离点则存入outline
    16.      if abs(x(i)-average2) >= standard2*3
    17.        outlier = [outlier x(i)];
    18.      end
    19. end
    20.

    ②替换异常值

    替换方法跟缺失值的替换一致，我们可以直接将其看作缺失值进行处理，替换方法如上。

    3、数据变换

    ①0-1标准化

    1. import numpy as np
    2. import pandas as pd
    3. import matplotlib.pyplot as plt
    4. # matplotlib inline
    5. # 数据标准化
    6. # （1）0-1标准化
    7. df = pd.DataFrame({"value1":np.random.rand(10)*20,
    8.                   'value2':np.random.rand(10)*100})
    9. print(df.head())
    10. print('------')
    11. # 创建数据
    12.  
    13. def data_norm(df,*cols):
    14.     df_n = df.copy()
    15.     for col in cols:
    16.         ma = df_n[col].max()
    17.         mi = df_n[col].min()
    18.         df_n[col + '_n'] = (df_n[col] - mi) / (ma - mi)
    19.     return(df_n)
    20. # 创建函数，标准化数据
    21. df_n = data_norm(df, 'value1',  'value2')
    22. print(df_n.head())#标准化数据

    ②z-score标准化

    1. import numpy as np
    2. import pandas as pd
    3. import matplotlib.pyplot as plt
    4.  
    5. df = pd.DataFrame({"value1": np.random.rand(10) * 100,
    6.                    'value2': np.random.rand(10) * 100})
    7. print(df.head())
    8. print('------')
    9.  
    10.  
    11. # 创建数据
    12.  
    13. def data_Znorm(df, *cols):
    14.     df_n = df.copy()
    15.     for col in cols:
    16.         u = df_n[col].mean()
    17.         std = df_n[col].std()
    18.         df_n[col + '_Zn'] = (df_n[col] - u) / std  # 平均值/标准差
    19.     return (df_n)
    20.  
    21.  
    22. # 创建函数，标准化数据
    23.  
    24. df_z = data_Znorm(df, 'value1', 'value2')
    25. u_z = df_z['value1_Zn'].mean()
    26. std_z = df_z['value1_Zn'].std()
    27. print(df_z)
    28. print('标准化后value1的均值为:%.2f, 标准差为：%.2f' % (u_z, std_z))
    29. # 标准化数据
    30. # 经过处理的数据符合标准正态分布，即均值为0，标准差为1
    31.  
    32. # 什么情况用Z-score标准化：
    33. # 在分类、聚类算法中，需要使用距离来度量相似性的时候，Z-score表现更好

    ③标准化的应用

    1. # 八类产品的两个指标value1，value2，其中value1权重为0.6，value2权重为0.4
    2. # 通过0-1标准化，判断哪个产品综合指标状况最好
    3. import numpy as np
    4. import pandas as pd
    5. import matplotlib.pyplot as plt
    6.  
    7. df = pd.DataFrame({"value1": np.random.rand(10) * 30,
    8.                    'value2': np.random.rand(10) * 100},
    9.                   index=list('ABCDEFGHIJ'))
    10.  
    11.  
    12. # print(df.head())
    13. # print('------')
    14. # 创建数据"
    15. def data_norm(df, *cols):
    16.     df_n = df.copy()
    17.     for col in cols:
    18.         ma = df_n[col].max()
    19.         mi = df_n[col].min()
    20.         df_n[col + '_n'] = (df_n[col] - mi) / (ma - mi)
    21.     return df_n
    22.  
    23.  
    24. df_n1 = data_norm(df, 'value1', 'value2')
    25. # 进行标准化处理
    26.  
    27. df_n1['f'] = df_n1['value1_n'] * 0.6 + df_n1['value2_n'] * 0.4
    28. df_n1.sort_values(by='f', inplace=True, ascending=False)
    29. df_n1['f'].plot(kind='line', style='--.k', alpha=0.8, grid=True)
    30. print(df_n1)
    31. # 查看综合指标状况

    4、数据离散化

    ①等宽法

    1. import numpy as np
    2. import pandas as pd
    3. import matplotlib.pyplot as plt
    4.  
    5. # 等宽法 → cut方法
    6. ages = [20, 22, 25, 27, 21, 23, 37, 31, 61, 45, 41, 32]
    7. # 有一组人员年龄数据，希望将这些数据划分为“18到25”,“26到35”,“36到60”,“60以上”几个面元，分成4个区间。
    8.  
    9. bins = [18, 25, 35, 60, 100]
    10. cats = pd.cut(ages, bins)
    11. print(cats)
    12. print(type(cats))
    13. print('____')
    14.  
    15. print(cats.codes, type(cats.codes))  # 0-3对应分组后的四个区间，用代号来注释数据对应区间，结果为ndarray；可以查看里边的等级
    16. print(cats.categories, type(cats.categories))  # 四个区间，结果为index
    17. print(pd.value_counts(cats))  # 按照区间计数
    18. print('-------')
    19. # cut结果含有一个表示不同分类名称的层级数组以及一个年龄数据进行标号的代号属性
    20.  
    21. print(pd.cut(ages, [18, 26, 36, 61, 100], right=False))
    22. print('-------')
    23. # 通过right函数修改闭端，默认为True
    24.  
    25. group_names = ['Youth', 'YoungAdult', 'MiddleAged', 'Senior']
    26. print(pd.cut(ages, bins, labels=group_names))
    27. print('-------')
    28. # 可以设置自己的区间名称，用labels参数
    29.  
    30.  
    31. df = pd.DataFrame({'ages': ages})
    32. group_names = ['Youth', 'YoungAdult', 'MiddleAged', 'Senior']
    33. s = pd.cut(df['ages'], bins)  # 也可以 pd.cut(df['ages'],5),将数据等分为5份
    34. df['label'] = s
    35. cut_counts = s.value_counts(sort=False)
    36. print(df)
    37. print(cut_counts)
    38. # 对一个Dataframe数据进行离散化，并计算各个区间的数据计数
    39.  
    40. plt.scatter(df.index, df['ages'], cmap='Reds', c=cats.codes)
    41. plt.grid()
    42. # 用散点图表示，其中颜色按照codes分类
    43. # 注意codes是来自于Categorical对象

    ②等频法

    1. # 等频法 → qcut方法
    2. import numpy as np
    3. import pandas as pd
    4. import matplotlib.pyplot as plt
    5. data = np.random.randn(1000)
    6. s = pd.Series(data)
    7. cats = pd.qcut(s,4)  # 按四分位数进行切割，可以试试 pd.qcut(data,10)
    8. print(cats.head())
    9. print(pd.value_counts(cats))
    10. print('------')
    11. # qcut → 根据样本分位数对数据进行面元划分，得到大小基本相等的面元，但并不能保证每个面元含有相同数据个数
    12. # 也可以设置自定义的分位数（0到1之间的数值，包含端点） → pd.qcut(data1,[0,0.1,0.5,0.9,1])
    13.  
    14. plt.scatter(s.index,s,cmap = 'Greens',c = pd.qcut(data,4).codes)
    15. plt.xlim([0,1000])
    16. plt.grid()
    17. # 用散点图表示，其中颜色按照codes分类
    18. # 注意codes是来自于Categorical对象

    三、数据降维

    1、主成分分析

    ①主成分分析简介

    主成分分析是一种降维算法，它能将多个指标转换为少数几个主成分，这些主成分是原始变量的线性组合，且彼此之间互不相关，其能反映出原始数据的大部分信息。一般来说，当研究的问题涉及到多变量且变量之间存在很强的相关性时，我们可考虑使用主成分分析的方法来对数据进行简化。

    ②主成分分析计算步骤

    标准化处理

    计算标准化样本的协方差矩阵

    计算R的特征值和特征向量

    计算主成分贡献率以及累计贡献率

    写出主成分

    根据系数分析主成分代表的意义

    ③案例分析

    主成分分析指标解释案例
    主成分分析的一大难点是指标意义模糊，难以解释，下面这个例子可以辅助理解。

    在这里插入图片描述

    在这里插入图片描述
    上表的累计贡献率 = 当前项贡献率 + 之前的累计贡献率。当累计贡献率 > 80%时，剩下的特征向量可以舍弃。

    在这里插入图片描述

    上面的分析需要一定的语言组织能力，也需要一定运气成分，若难以解释，或者强行解释，或者换方法。

    案例参考文章：原文链接：https://blog.csdn.net/qq1198768105/article/details/119898545

    在这里插入图片描述

    在这里插入图片描述

    在这里插入图片描述

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

     2、因子分析（FA）

    ①应用场景

    减少分析变量个数
    通过对变量间相关关系的探测，将原始变量分组，即将相关性高的变量分为一组，用共性因子来代替该变量
    使问题背后的业务因素的意义更加清晰呈现

    ②步骤

    1、选择分析的变量

    2、计算所选原始变量的相关系数矩阵

    3、提取公共因子

    4、因子旋转

    5、计算因子得分

    ③模型分析

     以废气排放为例说明主成分降维处理过程

    1、评价标准体系的构建

    2、因子分析

    因子分析首先将原始数据标准化处理,建立相关系数矩 阵并计算其特征值和特征向量,接着从中选择特征值大于等 于1的特征值个数为公共因子数,或者根据因子对X的累计贡献 率大于80%来确定公共因子,求得因子载荷矩阵, 后计算公因子得分和综合得分。

     观察相关系数矩阵表 2，可以发现所选取指标之间存在着一定的相关关系，其中 X3 和 X6、X4 和
    X7、 X7 和 X8 分别存在着较强的相关性，相关系数分别为 0.96、 0.96、0.91，这进一步验证了对所选指标做因子分析的科学性和必要性。计算相关系数矩阵的特征值、贡献率及累计贡献率如
    所示：

     确定提取的主成分个数可综合考虑 3 个方面:

    (1)提取的所有特征值大于某一特定特征值，一般特定值设为 1，本文同样以 1 为标准；

    (2) 提取的主成分的累计贡献率要大于 85%，即所提取的主成分要能够概括原有指标的绝大部分信息; 由表 3 可知，前 3 个主成分的累计贡献率已经达到了94.01%，满足按照累计贡献率大于 85% 确定主成分个数的原则；

    (3) 以做主成分分析时生成的碎石图 (Scree Plot) 做参考，碎石图是按照特征值大小排列的，以特征值为纵坐标、因子数为横坐标生成的主成分散点图，有明显的拐点，一般取拐点前所有的因子及拐点后第一个因子作为主成分 。
    观察图 4 可得，第一、第二个主成分的特征值较大，其余几个均较小，碎石图在第三个特征值出现
    拐点。
    根据上述分析，在本研究中选取前 3 个主成分对河南省水资源使用情况进行动态分析，从表 3
    和 4，我们可以看出前 3 个主成分已经能够概括绝大部分的原始信息，因此提取 3 个主成分因子是
    合理的. 提取 3 个主成分用于概括原有 10 个指标的绝大部分信息，这既达到了降维、简化的目的，又在一定程度上保证了后续研究结果能准确有效地反映出河南省水资源使用情况动态变化的基本特征。计算主成分的载荷矩阵, 主成分载荷是指提取的 3 大主成分与各变量指标之间的相关系数如表4：

     ④代码实现

    1. clc,clear
    2. r=[1.000 0.577 0.509 0.387 0.462
    3. 0.577 1.000 0.599 0.389 0.322
    4. 0.509 0.599 1.000 0.436 0.426
    5. 0.387 0.389 0.436 1.000 0.523
    6. 0.462 0.322 0.426 0.523 1.000];
    7. %下面利用相关系数矩阵求主成分解，val的列为r的特征向量，即主成分的系数
    8. [vec,val,con]=pcacov(r);%val为r的特征值，con为各个主成分的贡献率
    9. f1=repmat(sign(sum(vec)),size(vec,1),1); %构造与vec同维数的元素为±1的矩阵
    10. vec=vec.*f1; %修改特征向量的正负号，每个特征向量乘以所有分量和的符号函数值
    11. f2=repmat(sqrt(val)',size(vec,1),1);
    12. a=vec.*f2 %构造全部因子的载荷矩阵
    13. a1=a(:,1) %提出一个因子的载荷矩阵
    14. tcha1=diag(r-a1*a1') %计算一个因子的特殊方差
    15. a2=a(:,[1,2]) %提出两个因子的载荷矩阵
    16. tcha2=diag(r-a2*a2') %计算两个因子的特殊方差
    17. ccha2=r-a2*a2'-diag(tcha2) %求两个因子时的残差矩阵
    18. gong=cumsum(con) %求累积贡献率
    19.  
    20.  
    21.  
    22. clc,clear
    23. load data.txt; %把原始数据保存在纯文本文件data.txt中
    24. n=size(data,1);
    25. x=data(:,1:4); y=data(:,5); %分别提出自变量x和因变量y的值
    26. ——————————————————————————————————
    27. 如果不需要检验，则不需要把y列入原始数据中，把矩阵x的大小改变一下，以及下文中的m，m为原始数据中变量的个数。
    28. ——————————————————————————————————
    29. m=4;%m为变量的个数
    30. x=zscore(x); %数据标准化
    31. r=cov(x); %求标准化数据的协方差阵，即求相关系数矩阵
    32. [vec,val,con]=pcacov(r); %进行主成分分析的相关计算
    33. c=cumsum(con);
    34. i=1;
    35. while ((c(i)<90)&(con(i+1)>10))
    36.     i=i+1;
    37. end
    38. num=i;
    39. f1=repmat(sign(sum(vec)),size(vec,1),1);
    40. vec=vec.*f1; %特征向量正负号转换
    41. f2=repmat(sqrt(val)',size(vec,1),1);
    42. a=vec.*f2; %求初等载荷矩阵
    43. am=a(:,1:num); %提出num个主因子的载荷矩阵
    44. [b,t]=rotatefactors(am,'method', 'varimax'); %旋转变换,b为旋转后的载荷阵
    45. bt=[b,a(:,num+1:end)]; %旋转后全部因子的载荷矩阵
    46. contr=sum(bt.^2); %计算因子贡献
    47. rate=contr(1:num)/sum(contr); %计算因子贡献率
    48. fprintf('综合因子得分公式：F=');
    49. for i=1:num
    50.     fprintf('+%f*F%d',rate(i),i);
    51. end
    52. fprintf('\n');
    53. coef=inv(r)*b; %计算得分函数的系数
    54. coef=coef';
    55. for i=1:num
    56.     fprintf('各个因子得分函数为F%d=',i);
    57.     for j=1:m
    58.         fprintf('+(%f)*x_%d',coef(i,j),j);
    59.     end
    60.     fprintf('\n');
    61. end
    62. %如果仅仅因子分析，程序到此为止
    63. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    64. score=x*coef';%计算各个因子的得分
    65. weight=rate/sum(rate); %计算得分的权重
    66. Tscore=score*weight'; %对各因子的得分进行加权求和，即求各企业综合得分
    67. [STscore,ind]=sort(Tscore,'descend'); %对企业进行排序
    68. display=[score(ind,:)';STscore';ind']; %显示排序结果
    69. fprintf('排序结果如下：');
    70. for i=1:num
    71.     fprintf('第%d行为F%d得分,',i,i);
    72. end
    73. fprintf('第%d行为综合因子得分，第%d为原序列\n',num+1,num+2);
    74. disp(display);
    75. [ccoef,p]=corrcoef([Tscore,y]); %计算F与资产负债的相关系数
    76. [d,dt,e,et,stats]=regress(Tscore,[ones(n,1),y]);%计算F与资产负债的方程
    77. fprintf('因子分析法的回归方程为：F=%f+(%f*y)',d(1),d(2));
    78. if (stats(3)<0.05)%判断是否通过显著性检验的结果
    79.     fprintf('\n在显著性水平0.05的情况下，通过了假设检验。\n');
    80. else
    81.     fprintf('\n在显著性水平0.05的情况下，通不过假设检验。\n');
    82. end该MATLAB源代码的displsy为最终排序结果。

  • 相关阅读:
    ARM编程环境搭建教程
    Linux系统运维排故思路参考手册
    floyd 算法模板详解(适合新手)
    Vue3 快速入门和模板语法
    第一章：概述
    Hadoop笔记
    SSH项目集成DM达梦数据库
    Abnova ABCB10（人）重组蛋白说明书
    Pytorch实用教程：nn.LSTM内部是如何实现的
    【数据库系统概论】第二章关系数据库
  • 原文地址：https://blog.csdn.net/m0_58585940/article/details/127664157