古玻璃制品的成分分析与鉴别详解【国一，附完整代码】

声明：2024年数模国赛即将来临，为助力国赛和钉钉杯，我将重温22年小样本国赛C题和23年大样本国赛C题，给出详细思路和完整代码，供广大数模爱好者阅览，如需比赛指导，请联系文章底部卡片咨询。
未经允许，不得转载。——CSDN:川川菜鸟

文章目录

题目：22年古玻璃制品的成分分析与鉴别
- 完整代码+数据文件
数据检验与预处理
- 先验知识：成分数据
- 1、缺失值
- 2、异常值
- 3、成分数据
问题一：探究规律+统计量预测
- 1、表单 1 中表面风化与玻璃类型、颜色和纹饰的关系
- 2、有无风化化学成分含量的统计规律
- 3、预测风化前各化学成分含量
问题二：分类+敏感性分析
- 1、两种玻璃的分类规律
- 2、玻璃亚类划分
- 3、亚类划分结果的敏感性探究
问题三：预测表单3的玻璃类别
问题四：不同化学成分间的关联关系
总结

题目：22年古玻璃制品的成分分析与鉴别

丝绸之路是古代中西方文化交流的通道，其中玻璃是早期贸易往来的宝贵物证。现有一批我国古代玻璃制品的相关数据，考古工作者依据这些文物样品的化学成分和其他检测手段已将其分为高钾玻璃和铅钡玻璃两种类型。
在这里插入图片描述
　　附件表单 1 给出了这些文物的分类信息，附件表单 2 给出了相应的主要成分所占比例（空白处表示未检测到该成分）。这些数据的特点是成分性，即各成分比例的累加和应为 100%，但因检测手段等原因可能导致其成分比例的累加和非 100%的情况。本题中将成分比例累加和85%~105%之间的数据视为有效数据。基于上述背景及附件中数据，本文需要建立数学模型解决以下问题：
问题1
　　分析表单 1 中表面风化与玻璃类型、颜色和纹饰的关系；结合玻璃类型，分析表面有无风化与化学成分含量的关系，并依据风化点监测数据，预测风化前各化学成分含量。
问题2
　　分析高钾玻璃、铅钡玻璃的分类规律；对每个类别选择合适的化学成分进行亚类划分，并分析分类结果的敏感性和合理性。
问题3
　　分析表单 3 中未知类别玻璃文物的化学成分，鉴定类别，并分析分类结果的敏感性。
问题4
　　对不同类别玻璃文物，分析其化学成分间的关联关系，并比较不同类别间化学成分关联关系的差异性。
（上述是简化后的题目，本文在想做题的同学直接在官网全国大学生数学建模竞赛官网下载题目与附件数据）

完整代码+数据文件

这里我把文件和代码全部放出来，读者自己去下载运行测试，后续将主要对代码做一些分析，帮助大家理解。

链接：https://pan.baidu.com/s/1t1RGW7yGjMs7SSFOR97ypQ 
提取码：ydee

数据检验与预处理

在数据挖掘类题目中，做前期的数据处理是必要的。既为避免“分母出现0值”从而出现inf，这种离谱的计算结果做准备，也为纠正某些数据的天然分布，对后续问题的处理造成误差做准备。比如此次题目的数据：成分数据！！！

先验知识：成分数据

“成分数据”通常指的是描述一个对象或系统中各个成分的数据集合。这类数据广泛应用于化学、生物学、食品科学、环境科学等领域。成分数据具备的重要特征：
多变量性：成分数据通常包括多个变量，每个变量代表不同的成分。这些变量之间可能存在相关性，因为它们共同组成了一个整体。（对应题目附件Excel中表单2和表单3，每列化学成分就是一个变量）
约束性：成分数据的一个重要特征是它们受到约束。例如，在化学成分分析中，所有成分的总和通常约束为100%（在百分比表示时）。这种约束导致数据间存在内在的依赖关系。（题目要求成分数据总和在85%~105%，这也是我们判别数据是否为异常值的标准）
比例性：成分数据通常以比例或百分比表示，反映各成分在总体中的相对量。这要求分析方法能适应比例数据的特点，比如使用适合处理比例数据的统计技术。（如对数比变换（CLR、ALR、ILR）、Dirichlet回归、Beta回归等）
正性：所有成分数据都是非负的，因为它们表示的是量或比例。（这说明如果算出来的化学成分为负值，那么得到的答案是不可信的）
高度的协同变化：在某些情况下，如果一个成分的比例增加，其他成分的比例必然减少（数据处理和分析时可以采取特殊的方法，例如主成分分析（PCA）的变体——对应分析（CA））

1、缺失值

观察附件数据，表单1、2、3中均有空值。考虑到数据处理中分母为0与成分数据的特征，表单1缺失的颜色在用众数“浅蓝”填充，表单2的空值用最小比例0.04填充，0值用0.04填充。表单3中的空值也用最小比例0.04填充，0值用0.04填充。（在数据处理中，空值和0值是不一样的！）

import pandas as pd
df1 = pd.read_excel('附件.xlsx', '表单1')
column_data = df1.iloc[:, 3]  #第几列,注意索引从0开始为第一列
mode_value = column_data.mode()[0]  #获取并选择第一个众数
df1.iloc[:, 3] = column_data.fillna(mode_value) #使用众数填充
df2 = pd.read_excel('附件.xlsx', '表单2')
column_data = df2.iloc[:, :] 
df2.iloc[:, :] = column_data.fillna(0.04) #填充为最小比例0.04
df2.replace(0, 0.04, inplace=True)  #将0值替换为0.04
df3 = pd.read_excel('附件.xlsx', '表单3')
column_data = df3.iloc[:, :] 
df3.iloc[:, :] = column_data.fillna(0.04) #填充为最小比例0.04
df3.replace(0, 0.04, inplace=True) #将0值替换为0.04
with pd.ExcelWriter('附件_处理后.xlsx', engine='openpyxl') as writer:
    df1.to_excel(writer, sheet_name='表单1', index=False)
    df2.to_excel(writer, sheet_name='表单2', index=False)
    df3.to_excel(writer, sheet_name='表单3', index=False)

运行截图如下：
在这里插入图片描述

2、异常值

依据题意，我们将表单2中成分数据加和在<85%或>105%外的样本进行剔除，具体代码如下:

import pandas as pd
xls = pd.ExcelFile('附件_处理后.xlsx')
df1 = xls.parse('表单2')
columns_to_sum = df1.iloc[:, 1:15]  # 选取第2列到第15列进行求和
row_sums = columns_to_sum.sum(axis=1)
df1['Sum'] = row_sums
filtered_df = df1[(df1['Sum'] >= 85) & (df1['Sum'] <= 105)]  #筛选出正确的值
excluded_df = df1[~((df1['Sum'] >= 85) & (df1['Sum'] <= 105))] # 筛选出异常值并查看
print(excluded_df)  #打印被剔除的样本
with pd.ExcelWriter('附件_处理后.xlsx', mode='a', engine='openpyxl', if_sheet_exists="replace") as writer: 
    # 将筛选后的数据写入'表单2'工作表，如果工作表存在则替换
    filtered_df.to_excel(writer, sheet_name='表单2', index=False)  #保存在原来的表单中
    excluded_df.to_excel(writer, sheet_name='被剔除的样本', index=False)

筛选出的异常值如下，整理好的数据保存在“附件_处理后.xlsx”中。
在这里插入图片描述

3、成分数据

由于成分数据的全部信息蕴含在各成分比例的相对大小中，仅考察某成分比例的绝对大小具有误导性，因此本文采用对数比变换CLR对附件表单2中，呈现右偏分布的数据进行相应转换。读者也可采用ALR、ILR进行转换。这里有一个注意点：通常做变换是对特征即一列进行转换，但对于成分数据，是对样本即一行进行数据变换。

#中心对数处理
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
def clr(dataMatrix):
    geometricMean = np.exp(np.mean(np.log(dataMatrix), axis=1))
    clrTransformedData = np.log(dataMatrix / geometricMean[:, np.newaxis])
    return clrTransformedData
df = pd.read_excel('对数比变换流程.xlsx')
df1 = df.iloc[:, 0:14]
dataMatrix = df1.values
clrTransformedData = clr(dataMatrix)
clrTransformedTable = pd.DataFrame(clrTransformedData, columns=df1.columns)
clrTransformedTable.to_excel('对数变换比流程结束.xlsx', index=False)
print("中心对数比变换处理完毕")

代码运行结果如下：
在这里插入图片描述
这里绘制中心对数比变换前后的核密度图，观察转换效果。代码如下：

import pandas as pd
import matplotlib.pyplot as plt
df = pd.read_excel('附件_处理后 _转换后.xlsx', '表单2')
df1 = df.iloc[:, 1:15]
df2 = df1.iloc[28]        #探究不同采样点不断更换行数即可
plt.style.use('seaborn-poster')  #seaborn-poster、#classic #不同风格
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']  #解决中文问题
plt.rcParams['axes.unicode_minus'] = False  #解决负号问题
df2.plot(kind='hist', bins=5, color='steelblue', edgecolor='black', density=True)
df2.plot(kind='kde', color='red')
plt.xlabel('化学成分含量比例')
plt.ylabel('密度')
plt.title('26采样点转换后化学成分分布')
plt.show()

绘制完成的核密度图如下：

07号文物转换前

18号文物转换前

26文物转换前

07号文物转换后

18号文物转换后

26号文物转换后

问题一：探究规律+统计量预测

数据源上，问题一涉及到附件中表单1和表单2。第一小问分析我们分析表面风化与玻璃类型、颜色及纹饰的关系。第二小问要求我们结合玻璃类型，分析表面有无风化与化学成分含量的关系，并预测风化前各化学成分含量。

1、表单 1 中表面风化与玻璃类型、颜色和纹饰的关系

首先观察表单1中数据，如下图所示。除去“编号”列，纹饰、类型、颜色、表面风化四列均为定类数据，因此我们要考虑的，是分析定类数据之间关系的的方法。比如卡方独立检验、Cramer的V统计量、Fisher精确检验、相关系数、Logistic回归。在这里插入图片描述
这是我们采用excel的数据透视表对表单1中数据进行特征观察，汇总列联表如下图所示：
请添加图片描述

可以发现，样本总数均在50以上，但其中部分数据点，存在0点。根据方法的特征，对于样本总量>40，列联表中数据均>5的类型与表面是否风化的关系，采用卡方独立检验。对于列联表中部分数据<5的纹饰、颜色与表面是否风化的关系，采用Fisher-Freeman-Halton检验。

下面是关于卡方检验和Fisher精确检验的两段介绍：

卡方独立检验：用于评估两个定类（名义或序数）变量之间的独立性。如果检验结果显示变量不独立，那么这意味着这两个变量之间可能存在某种关系或依赖。
Fisher精确检验：对于样本量较小或数据极端不平衡的列联表，Fisher精确检验提供了一个替代方案。这种方法不依赖于大样本理论，因此在小样本情况下比卡方检验更为精确。（这里使用改进的Fisher精确检验）

基本实施均是构建列联表，然后根据相应方法内容计算，这里我给出数据透视图整理出相应代码及结果：

#卡方检验代码
import numpy as np
from scipy.stats import chi2_contingency
observed = np.array([[6, 12],   #频数表
                     [28, 12],
                      ])
chi2, p, dof, expected = chi2_contingency(observed)
print(f"卡方统计量: {chi2}")
print(f"p值: {p}")
print(f"自由度: {dof}")
print(f"期望频数表:\n{expected}")
alpha = 0.05  #95%的置信度
if p < alpha:
    print("拒绝原假设，行和列之间存在显著关联。")
else:
    print("不能拒绝原假设，行和列之间没有显著关联。")

#执行Fisher-Freeman-Halton检验
import numpy as np
import statsmodels.api as sm
contingency_table = np.array([[11, 11],   #列联表
                              [6, 0],
                              [17, 13]
                      ])
table = sm.stats.Table(contingency_table)
fisher_result = table.test_nominal_association()
print(f"p值: {fisher_result.pvalue}")
alpha = 0.05  #百分之95的置信度
if fisher_result.pvalue < alpha:
    print("拒绝原假设，行和列之间存在显著关联。")
else:
    print("不能拒绝原假设，纹饰与表面风化没有显著关联。")

置信度均取95%，结果如下图所示：
在这里插入图片描述

2、有无风化化学成分含量的统计规律

数据预处理，题目要求结合玻璃类型探究有无风化化学成分含量的统计规律，这就要求我们将附件数据中的表单1和表单2的数据整合，综合探究分析。首先使用excel处理表单2”文物采样点“列数据，只剩下数字编号并修改列名为“文物编号”复制到新表table1.xlsx中，再在表单1中挑出”文物编号“和“表面风化”两列保存在table2.xlsx中。下述是将table1.xlsx和table2.xlsx 整合的代码。

#将两张表的对应关系匹配
import pandas as pd
table1 = pd.read_excel('table2.xlsx',)
table2 = pd.read_excel('table1.xlsx')
result = pd.merge(table1, table2, on='文物编号', how='left')
result.to_excel('table1_updated.xlsx', index=False)
print("匹配完成，结果已保存到table1_updated.xlsx")

之后依据解释和上述代码，对“文物编号”和“玻璃类型”列做一样的操作，将两列结果复制到“附件_处理后 _转换后.xlsx”末尾两列，最终得到的结果如下所示：
请添加图片描述
下面依据“类型”列，将新整合的表单2分割成“铅钡.xlsx”和“高钾.xlsx”，再在“高钾.xlsx”和““高钾.xlsx”数据集上再次运行代码，得到**“高钾_风化.xlsx”、“高钾_无风化.xlsx”、“铅钡_风化.xlsx”、“铅钡_无风化.xlsx”**依次进行探究。代码如下：

#依据列的类别划分数据集
import pandas as pd
file_path = '附件_处理后 _转换后.xlsx'
df = pd.read_excel(file_path,'表单2')
column_to_split = '类型'
unique_values = df[column_to_split].unique()
for value in unique_values:
    df_subset = df[df[column_to_split] == value]
    new_file_path = f'{value}.xlsx'
    df_subset.to_excel(new_file_path, index=False)
    print(f'已保存: {new_file_path}')
print('所有文件已保存完成。')

最终得到结果如下所示：
请添加图片描述
官方还给了大家一点小坑，表单1中风化的文物在表单2中采样点处是无风化的，如下图所示，这就需要大家细致的阅读题目，以及对数据敏锐的洞察能力。数据较少，人工修改之后继续进行下一步操作。

请添加图片描述
为考虑玻璃类型，探究有无风化化学成分含量的统计规律，我们绘制箱型图，观测两种玻璃关于某种化学成分风化前后含量的差异，代码如下：

#绘制箱型图（多子图）的代码
import pandas as pd
import matplotlib.pyplot as plt
files = ['高钾_风化.xlsx', '高钾_无风化.xlsx', '铅钡_风化.xlsx', '铅钡_无风化.xlsx']
labels = ['高钾_风化', '高钾_无风化', '铅钡_风化', '铅钡_无风化']
dfs = [pd.read_excel(file).iloc[:, 1:15] for file in files]
column_names = dfs[0].columns
plt.style.use('seaborn-darkgrid')
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei'] # 处理中文乱码问题
plt.rcParams['axes.unicode_minus'] = False # 处理坐标轴负号问题
fig, axes = plt.subplots(nrows=7, ncols=2, figsize=(15, 30))
axes = axes.flatten()
for i in range(dfs[0].shape[1]): # 假设每个文件的第2到第15列都有相同的化学成分
    data = [df.iloc[:, i] for df in dfs]
    ax = axes[i]
    ax.boxplot(data, labels=labels)
    ax.set_title(f'{column_names[i]}', fontsize=10)
    ax.set_xlabel('样本分类', fontsize=10)
    ax.set_ylabel('含量', fontsize=10)

#ax.tick_params(axis='x', rotation=45)
for j in range(i + 1, len(axes)):
    fig.delaxes(axes[j])
plt.tight_layout()
plt.subplots_adjust(top=0.96) # 调整顶部距离以防止标题重叠
plt.show()

绘制图片如下所示：
请添加图片描述

3、预测风化前各化学成分含量

由于没有详细的每种玻璃的每个文物风化前后的化学成分差异，这里我们使用简单的均值统计量进行风化前化学成分含量的预测。
这里使用excel进行简单的数据处理，直接计算“高钾_风化.xlsx”、“高钾_无风化.xlsx”、“铅钡_风化.xlsx”、“铅钡_无风化.xlsx”中每列化学成分的均值及均值间的差。最终如下所示：
将上述差值在风化的数据上直接相加就可以了。

问题二：分类+敏感性分析

分析高钾玻璃、铅钡玻璃的分类规律要求我们根据化学成分对玻璃分类；还要选择合适的化学成分划分亚类，涉及到特征选择。具体的划分方法是一个表达式或者一种方法，及划分结果，并分析分类结果的合理性和敏感性。

1、两种玻璃的分类规律

为快速判别使用哪种分类模型，依据化学成分对铅钡玻璃的类型进行高效分类。这里借助matlab的分类工具箱，进行基本分类。

首先声明，为避免风化前后化学成分含量的差异对类别判定造成干扰，我们采用中心对数比变换后的无风化的玻璃成分数据集。关于数据集的得出重新跑一遍问题一第二小问分割数据集的代码。导入数据为经过CLR变换后的化学成分数据。具体步骤如下所示：
①切换工作路径，黏贴即可
②点击导入数据
请添加图片描述
③点击导入的数据
④打开导入的数据

⑤选中要导入的列，命名为S1，并点击导入
在这里插入图片描述
⑥打开APP，点击分类学习器

⑦导入变量，开始会话。

⑧点击全部，训练所有模型

⑨最终得到的分类准确率如下图所示：
请添加图片描述
可以发现，在所有模型中，matlab给出了准确率最高的模型为线性核支持向量机，准确率达到97.2%，这说明中心对数比变换后的数据集具有线性可分的性质，恰好后文要求探究分类规律，线性核的核函数作为分类规律正合适。因此我决定以线性核支持向量机为基本模型，后续对此进行改进。这里给出支持向量机相应的特征图和分类正确率图：

分类混淆矩阵

ROC曲线图

以下是matlab导出的线性核函数，作为铅钡玻璃和高钾玻璃的分类规律：

    y = -0.3993 + (0.0871 * SiO2) + (-0.0041 * Na2O) + (1.1271 * K2O) + (0.5476 * CaO) + (0.4565 * MgO) + (0.3753 * Al2O3) + (0.4789 * Fe2O3) + (0.7290 * CuO) + (-1.4872 * PbO) + (-1.2129 * BaO) + (0.3107 * P2O5) + (-0.5440 * SrO) + (0.2056 * SnO2) + (0.0196 * SO2)

这里以线性核函数为分类界限，一侧是高钾玻璃，另一侧是铅钡玻璃。即大于0为高钾，小于0为铅钡。

2、玻璃亚类划分

这里所用的数据，是中心对数比变换后的无风化的两种类型玻璃，即前文得到的“铅钡_无风化.xlsx”，“高钾_无风化.xlsx”。这里我们首先进行层次聚类，绘制系谱图。代码如下所示：

#绘制层次聚类系谱图
import pandas as pd
import matplotlib.pyplot as plt
import scipy.cluster.hierarchy as sch
file_path = '铅钡_无风化.xlsx'
data = pd.read_excel(file_path)
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']  #解决中文问题
plt.rcParams['axes.unicode_minus'] = False  #解决负号问题

print(data.head())
#选择要聚类的行
data_for_clustering = data.iloc[:, 1:15]
print("Selected data shape:", data_for_clustering.shape)
# 展现聚类
linked = sch.linkage(data_for_clustering, method='ward')
plt.figure(figsize=(10, 7))  #图片大小
sch.dendrogram(linked,
               orientation='top',
               distance_sort='descending',
               show_leaf_counts=True)
plt.title('铅钡层次聚类系谱图')
plt.xlabel('样本')
plt.ylabel('距离')
plt.show()

高钾玻璃和铅钡玻璃的层次聚类系谱图如下所示：

高钾层次聚类系谱图

铅钡层次聚类系谱图

对于为无风化_高钾类，我们首先根据层次聚类划分出大体的类别，选择高度为11划分出3类，**标记为1类，2类，3类。**之后为找到具体的亚类划分规律，我们使用**ReliefF算法**筛选出重要特征来简化问题，承接上文分类学习器，操作步骤如下图所示：

请添加图片描述
如下图所示，最终得到Na2O,P2O5,BaO,SO2,PbO5个特征。

依据5个特征，采用matlab分类工具箱进行快速训练，最终得到如下结果。发现核逼近分类器——SVM核的效果最好，我们决定基于SVM核的核逼近分类器后续对高钾玻璃亚类进行划分。
请添加图片描述
同样的，我们也对高钾玻璃实行如此操作，得到的重要特征如下所示：

也选取前5个特征，采用matlab分类工具箱进行快速训练，最终得到如下结果。发现随机森林模型的效果最好，我们决定利用这里训练出的随机森林模型对后续对铅钡玻璃亚类进行划分。

请添加图片描述

3、亚类划分结果的敏感性探究

这里使用第一问还原后的成分数据结果（未经中心对数比变换），重新使用层次聚类模型、核逼近器——SVM核、随机森林模型，进行亚类划分。再使用第二小问训练出的核逼近器和随机森林模型，进行亚类别预测，得到二者结果相差不大。说明所建模型对包含噪声的成分还原后的数据不敏感，仍发挥功用，说明我们所建模型具有鲁棒性。

问题三：预测表单3的玻璃类别

为预测表单3的玻璃类别，首先确定目的，不仅预测A1-A6的大类别，是高钾玻璃还是铅钡玻璃，其次要预测A1-A6分别属于哪个亚类。

导出第二问第一小问的线性支持向量机模型，导出名称为trainedModel，如下图所示：
在这里插入图片描述
之后整理表单3的数据集，将风化玻璃的化学成分依据均值还原为无风化玻璃的化学成分，并导入MATLAB命名为S2，使用语法**[yfit,scores]=trainedModel.predictFcn(S2)** 来预测，得到结果如下所示：
请添加图片描述
可以看到，我们的预测结果和官方评阅要点一致，全部预测准确。在这里插入图片描述
为对分类结果的敏感性进行分析，我们使用化学成分未还原的数据集进行训练，最终分类结果和上图相同，说明结果对噪声不敏感，甚至对风化后的玻璃文物，仍能有效分类。

为进行相关亚类划分，依照得出的高钾、铅钡两种玻璃类型，我们根据两种类型还原风化前化学成分含量，之后将高钾玻璃分类的基于SVM核的核逼近分类器导出模型trainedModel3,将铅钡玻璃分类的基于随机森林导出模型trainedModel4，对表单3中的数据进行预测，操作流程及结果如下所示，其中S3g是表单3中高钾玻璃的数据行，S4是表单3中铅钡玻璃的数据行。最终得出的数字是所属亚类。
高钾玻璃亚类预测
铅钡玻璃亚类预测

问题四：不同化学成分间的关联关系

探究不同类别的玻璃文物的化学成分关联关系，需要对玻璃文物分类，直接对第2问得到的”铅钡.xlsx“和"高钾.xlsx"数据集进行处理。由于给出的数据只有一些化学成分，探究其内在联系，可以划入无监督问题。

思考基于奇异值分解的主成分分析，通过矩阵分解的方法，确切探究数据间关系，结合主成分分析和山崖落石图，可得到主要成分，再通过主成分因子载荷矩阵可得到主要成分的起主要贡献的化学成分。最后通过双坐标图，可探究风化前后或整体化学成分的关联关系。

首先绘制山崖落石图，确定每种类型玻璃的主成分选取几个合适，对应代码如下所示：

#绘制山崖落石图
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# 读取Excel文件
file_path = '高钾.xlsx'  # Excel文件的路径
sheet_name = '表单2'  # Excel工作表的名称
df = pd.read_excel(file_path)
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']  #解决中文问题
plt.rcParams['axes.unicode_minus'] = False  #解决负号问题
# 选择第2到第15行的数据，并确保数据是数值类型
data_matrix = df.iloc[:, 2:15]
# 对矩阵转置后进行奇异值分解
U, s, Vt = np.linalg.svd(data_matrix.T)
# 将奇异值转换为对角矩阵
Sigma = np.diag(s)
# 山崖落石图（Scree Plot）
plt.figure(figsize=(8, 5))
plt.plot(s**2 / np.sum(s**2), 'o-', linewidth=2, color='blue')
plt.title('高钾玻璃的山崖落石图')
plt.xlabel('Component Number')
plt.ylabel('Eigenvalue (Explained Variance Ratio)')
plt.grid(True)
# 显示图表
plt.show()

绘制山崖落石图如下所示，图像的拐点处即为合适的主成分个数。

请添加图片描述

可以看到两幅图中，在主成分个数为4的时候，可以解释绝大部分方差。因此我们均选择四个主成分进行后续操作。

其次绘制化学成分的主成分载荷矩阵图，代码如下所示：

#载荷矩阵图代码
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
# 读取Excel文件
file_path = '铅钡.xlsx'  # Excel文件的路径
sheet_name = '表单2'  # Excel工作表的名称
df = pd.read_excel(file_path)
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']  # 解决中文问题
plt.rcParams['axes.unicode_minus'] = False  # 解决负号问题
# 选择第2到第15行的数据，并确保数据是数值类型
data_matrix = df.iloc[:, 2:15].values
# 对矩阵转置后进行奇异值分解
U, s, Vt = np.linalg.svd(data_matrix.T)
# 取前四个主成分的载荷
loadings = U[:, :4]
# 将载荷矩阵转换为DataFrame并重塑为长格式
loadings_df = pd.DataFrame(loadings, columns=['PC1', 'PC2', 'PC3', 'PC4'])
loadings_long = loadings_df.melt(var_name='Principal Component', value_name='Loading', ignore_index=False)
loadings_long['Variable'] = np.tile(df.columns[2:15], 4)

# 绘制气泡图
plt.figure(figsize=(8, 8))
bubble_plot = sns.scatterplot(data=loadings_long, x='Principal Component', y='Variable', size='Loading', hue='Loading',
                              sizes=(600, 600), palette='coolwarm', legend=None, edgecolor='w', linewidth=0.5)
# 调整气泡的颜色条
norm = plt.Normalize(loadings_long['Loading'].min(), loadings_long['Loading'].max())
sm = plt.cm.ScalarMappable(cmap="coolwarm", norm=norm)
sm.set_array([])
bubble_plot.get_figure().colorbar(sm, orientation="vertical", fraction=0.046, pad=0.04)
# 设置标签和标题
bubble_plot.set_title('铅钡玻璃的载荷矩阵气泡图', fontsize=16)
bubble_plot.set_xlabel('主成分', fontsize=14)
bubble_plot.set_ylabel('变量', fontsize=14)
# 调整布局
plt.tight_layout()
#plt.grid(True)
# 显示图表
plt.show()

绘制主成分载荷矩阵的图片如下所示，可以看到铅钡玻璃的第一主成分起主要贡献的是PbO和BaO,第二主成分起主要贡献的是Na2O，第三主成分起主要贡献的是SO2和CuO,第四主成分起主要贡献的是MgO。
高钾玻璃的第一主成分起主要贡献的是SnO2和SrO，第二主成分起主要贡献的是BaO和PbO，第三主成分起主要贡献的是MgO和P2O5，第四主成分起主要贡献的是PbO。
请添加图片描述

最后绘制高钾玻璃和铅钡玻璃风化前后的双标图，和单个类型玻璃整体的双标图，图中的颜色阴影为风化前后的置信区域，观察化学成分间的关联关系。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from matplotlib.patches import Ellipse

# 读取Excel文件
file_path1 = '铅钡_无风化.xlsx'  # Excel文件的路径
file_path2 = '铅钡_风化.xlsx'  # Excel文件的路径
sheet_name = '表单2'  # Excel工作表的名称

df1 = pd.read_excel(file_path1)
df2 = pd.read_excel(file_path2)

plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']  # 解决中文问题
plt.rcParams['axes.unicode_minus'] = False  # 解决负号问题

# 选择第2到第15列的数据，并确保数据是数值类型
data_matrix1 = df1.iloc[:, 2:15].values
data_matrix2 = df2.iloc[:, 2:15].values

# 进行主成分分析
pca = PCA(n_components=2)
principal_components1 = pca.fit_transform(data_matrix1)
principal_components2 = pca.transform(data_matrix2)  # 使用同一个PCA对象转换第二个数据集
loadings = pca.components_.T

# 绘制双标图
plt.figure(figsize=(10, 8))

# 绘制样本点
plt.scatter(principal_components1[:, 0], principal_components1[:, 1], color='red', edgecolor='k', s=50, label='无风化',
            alpha=0.6)
plt.scatter(principal_components2[:, 0], principal_components2[:, 1], color='blue', edgecolor='k', s=50, label='风化',
            alpha=0.6)

for i, (x, y) in enumerate(principal_components1):
    plt.text(x, y, str(i + 1), color='red', fontsize=9, ha='right')
for i, (x, y) in enumerate(principal_components2):
    plt.text(x, y, str(i + 1 + len(principal_components1)), color='blue', fontsize=9, ha='right')

# 绘制变量载荷
for i, (x, y) in enumerate(loadings):
    plt.arrow(0, 0, x, y, color='green', alpha=0.8, head_width=0.05, linewidth=1.5)
    plt.text(x * 1.15, y * 1.15, df1.columns[i + 2], color='green', ha='center', va='center', fontsize=10)


# 绘制阴影区域（椭圆）
def plot_ellipse(data, ax, n_std=2.0, facecolor='none', **kwargs):
    #绘制一个包含给定数据点的置信椭圆。
    cov = np.cov(data, rowvar=False)
    mean = np.mean(data, axis=0)

    pearson = cov[0, 1] / np.sqrt(cov[0, 0] * cov[1, 1])
    ell_radius_x = np.sqrt(1 + pearson)
    ell_radius_y = np.sqrt(1 - pearson)
    ellipse = Ellipse((0, 0), width=ell_radius_x * 2, height=ell_radius_y * 2, facecolor=facecolor, **kwargs)

    scale_x = np.sqrt(cov[0, 0]) * n_std
    mean_x = mean[0]
    scale_y = np.sqrt(cov[1, 1]) * n_std
    mean_y = mean[1]

    transf = plt.matplotlib.transforms.Affine2D() \
        .rotate_deg(45) \
        .scale(scale_x, scale_y) \
        .translate(mean_x, mean_y)

    ellipse.set_transform(transf + ax.transData)
    return ax.add_patch(ellipse)


ax = plt.gca()
plot_ellipse(principal_components1, ax, n_std=2.0, edgecolor='red', facecolor='red', alpha=0.2)
plot_ellipse(principal_components2, ax, n_std=2.0, edgecolor='blue', facecolor='blue', alpha=0.2)

# 设置标签和标题
plt.xlabel(f'主成分 1 ({pca.explained_variance_ratio_[0]:.2%})', fontsize=14)
plt.ylabel(f'主成分 2 ({pca.explained_variance_ratio_[1]:.2%})', fontsize=14)
plt.title('铅钡玻璃的主成分双标图', fontsize=16)

# 设置坐标轴范围
plt.xlim([-0.95, 0.95])
plt.ylim([-0.9, 0.9])

# 添加网格和图例
plt.grid(True, linestyle='--', alpha=0.7)
plt.axhline(0, color='black', linewidth=0.5)
plt.axvline(0, color='black', linewidth=0.5)
plt.legend(loc='upper right', fontsize=12)

# 显示图表
plt.tight_layout()
plt.show()
"""
"""
#绘制主成分分析的双坐标图
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA

file_path = '高钾.xlsx'  # Excel文件的路径
sheet_name = '表单2'  # Excel工作表的名称
df = pd.read_excel(file_path)
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']  # 解决中文问题
plt.rcParams['axes.unicode_minus'] = False  # 解决负号问题

data_matrix = df.iloc[:, 2:15].values

pca = PCA(n_components=2)  #主成分分析
principal_components = pca.fit_transform(data_matrix)
loadings = pca.components_.T
plt.figure(figsize=(10, 8))  #主成分分析

plt.scatter(principal_components[:, 0], principal_components[:, 1], edgecolor='k', s=50)
for i, (x, y) in enumerate(principal_components):
    plt.text(x, y, str(i+1), color='k')
for i, (x, y) in enumerate(loadings):  #绘制变量载荷
    plt.arrow(0, 0, x, y, color='r', alpha=0.5, head_width=0.05)
    plt.text(x * 1.15, y * 1.15, df.columns[i+2], color='r', ha='center', va='center')

plt.xlabel(f'主成分 1 ({pca.explained_variance_ratio_[0]:.2%})')  #主成分
plt.ylabel(f'主成分 2 ({pca.explained_variance_ratio_[1]:.2%})')  #主成分
plt.title('高钾玻璃的主成分双标图')
plt.xlim([-0.8, 0.8])  #
plt.ylim([-0.8, 0.8])
plt.grid(True)   #添加网格和显示图表
plt.axhline(0, color='black', linewidth=0.5)
plt.axvline(0, color='black', linewidth=0.5)
plt.show()

高钾玻璃的双标图如下所示：
请添加图片描述

铅钡玻璃的双标图如下所示：
请添加图片描述

对于上述图片的分析，两个箭头间的连线短，箭头近，说明比例变量的对数比接近常数。如果连线很长，说明对数比是高度变化的。两个箭头间的角度相当于两个对数比间的相关系数 ，不相关的对数比的箭头正交，0度或180度箭头的对应对数比，完全线性相关。依据此原理分析即可。

总结

本文是对2022年国赛C题的全面复盘。首先对缺失值、异常值进行简单处理，再对成分数据进行中心对数比转换。之后采用卡方独立检验和Fisher-Freeman-Halton检验对定类数据进行处理，绘制箱型图探究风化前后化学成分的差异，再利用平均值统计量来还原风化前化学成分含量。
其次利用线性核的支持向量机进行玻璃大类划分。接着进行特征选择后，利用SVM核——核逼近分类器和随机森林进行高钾和铅钡玻璃的亚类划分。同时对人工还原的玻璃文物进行类别鉴定验证本文所建模型的鲁棒性。
接着重新运行问题二的所有模型，得到鉴定其所属类别和亚类，最终结果和官方评阅要点相同。
最后利用基于奇异值分解的主成分分析进行化学成分间关联关系探究。首先绘制山崖落石图选定合适的主成分个数，再绘制载荷矩阵气泡图观察每个主成分里起主要贡献的化学成分，最后绘制风化前后双标图和整体双标图，挖掘各化学成分间关联关系。