为了解决该问题，采用主成分分析（PCA）和多元线性回归模型来进行建模。PCA用于降维，帮助我们识别与玻璃类型、风化状态、纹饰、颜色相关的主要成分，进而简化数据结构。多元线性回归则用于基于风化后的成分预测风化前的化学成分含量。
主成分分析（PCA） (1) 数据预处理首先，将每个样本的化学成分比例矩阵定义为X∈Rn × pX \in \mathbb{R}^{n \times p}X∈Rn×p，其中nnn 是样本数量，ppp 是化学成分的数量。由于数据中的成分累加和并不总是100%，对数据进行归一化处理。假设第iii 个样本的化学成分比例为x i =[xi 1,xi 2,…,xi p]\mathbf{x}_i = [x_{i1}, x_{i2}, \dots, x_{ip}]xi=[xi1,xi2,…,xip]，则归一化后的成分表示为：
xi j′ = xij∑j=1pxijx'_{ij} = \frac{x_{ij}}{\sum_{j=1}^{p} x_{ij}}xij′=∑j=1pxijxij
确保成分比例的总和接近1。
(2) 协方差矩阵计算在完成数据归一化后，构建数据矩阵X ′ X'X′ 并计算其协方差矩阵Σ\SigmaΣ：
Σ=1n − 1∑i = 1n (x i ′ − x ˉ′ )(x i ′ − x ˉ′ ) ⊤ \Sigma = \frac{1}{n-1} \sum_{i=1}^{n} (\mathbf{x}'_i - \bar{\mathbf{x}}') (\mathbf{x}'_i - \bar{\mathbf{x}}')^\topΣ=n−11∑i=1n(xi′−xˉ′)(xi′−xˉ′)⊤
其中， x ˉ′ \bar{\mathbf{x}}'xˉ′ 是样本的均值向量。协方差矩阵Σ\SigmaΣ 用于捕捉成分之间的相关性。
(3) 特征值分解接下来，进行特征值分解（EVD），即：
Σ=VΛV ⊤ \Sigma = V \Lambda V^\topΣ=VΛV⊤
其中，Λ=diag(λ 1 ,λ 2 ,…,λ p )\Lambda = \text{diag}(\lambda_1, \lambda_2, \dots, \lambda_p)Λ=diag(λ1,λ2,…,λp) 为协方差矩阵的特征值，λ 1 ≥λ 2 ≥⋯≥λ p \lambda_1 \geq \lambda_2 \geq \dots \geq \lambda_pλ1≥λ2≥⋯≥λp，V=[v 1 ,v 2 ,…,v p ]V = [\mathbf{v}_1, \mathbf{v}_2, \dots, \mathbf{v}_p]V=[v1,v2,…,vp] 是协方差矩阵的特征向量。特征向量v i \mathbf{v}_ivi 对应主成分方向，而特征值λ i \lambda_iλi 则表明该主成分解释的方差大小。
(4) 主成分选择为了简化模型，只选择前kkk 个特征值较大的主成分。累积贡献率表示为：
累积贡献率=∑i=1kλi∑i=1pλi\text{累积贡献率} = \frac{\sum_{i=1}^{k} \lambda_i}{\sum_{i=1}^{p} \lambda_i}累积贡献率=∑i=1pλi∑i=1kλi
通常选择使累积贡献率达到80%-90%的kkk 值。选定的主成分矩阵为Z=XV k Z = XV_kZ=XVk，其中V k V_kVk 是前kkk 个主成分方向组成的矩阵。
(5) PCA结果分析将PCA降维后的数据ZZZ 与玻璃类型、颜色、纹饰进行关联分析。使用PCA降维后的二维或三维散点图观察不同玻璃类型（如高钾玻璃、铅钡玻璃）的分布，以及不同风化程度与颜色的关联性。
多元线性回归模型 (1) 模型设定我们要基于风化后的化学成分预测其风化前的成分。假设风化后的成分为x post =[xpost , 1,xpost , 2,…,xpost , p]\mathbf{x}_{\text{post}} = [x_{\text{post},1}, x_{\text{post},2}, \dots, x_{\text{post},p}]xpost=[xpost,1,xpost,2,…,xpost,p]，风化前的成分为x pre =[xpre , 1,xpre , 2,…,xpre , p]\mathbf{x}_{\text{pre}} = [x_{\text{pre},1}, x_{\text{pre},2}, \dots, x_{\text{pre},p}]xpre=[xpre,1,xpre,2,…,xpre,p]，则可以构建多元线性回归模型：
x pre =Bx post +ϵ\mathbf{x}_{\text{pre}} = \mathbf{B} \mathbf{x}_{\text{post}} + \mathbf{\epsilon}xpre=Bxpost+ϵ
其中，B∈Rp × p\mathbf{B} \in \mathbb{R}^{p \times p}B∈Rp×p 为回归系数矩阵，ϵ\mathbf{\epsilon}ϵ 为残差项，表示模型误差。
(2) 回归系数的估计为了估计回归系数矩阵B\mathbf{B}B，采用最小二乘法（OLS, Ordinary Least Squares），即：
B ^ =arg min B ∑i = 1n ∥xpre , i−Bxpost , i∥ 2 2 \hat{\mathbf{B}} = \arg \min_{\mathbf{B}} \sum_{i=1}^{n} \|\mathbf{x}_{\text{pre},i} - \mathbf{B} \mathbf{x}_{\text{post},i}\|_2^2B^=argminB∑i=1n∥xpre,i−Bxpost,i∥22
其解析解为：
B ^ =(X post ⊤ X post )− 1X post ⊤ X pre \hat{\mathbf{B}} = (\mathbf{X}_{\text{post}}^\top \mathbf{X}_{\text{post}})^{-1} \mathbf{X}_{\text{post}}^\top \mathbf{X}_{\text{pre}}B^=(Xpost⊤Xpost)−1Xpost⊤Xpre
其中，X post \mathbf{X}_{\text{post}}Xpost 和X pre \mathbf{X}_{\text{pre}}Xpre 分别为风化后和风化前成分的数据矩阵。
(3) 模型评估与预测通过交叉验证评估模型的预测能力。将模型应用于风化后的新数据，预测其风化前的化学成分：
x ^pre =B ^ x post \hat{\mathbf{x}}_{\text{pre}} = \hat{\mathbf{B}} \mathbf{x}_{\text{post}}x^pre=B^xpost
对比实际数据与预测结果，计算预测误差，例如使用均方误差（MSE）评估：
MSE=1 n ∑i = 1n ∥ x ^ pre , i−xpre , i∥ 2 2 \text{MSE} = \frac{1}{n} \sum_{i=1}^{n} \|\hat{\mathbf{x}}_{\text{pre},i} - \mathbf{x}_{\text{pre},i}\|_2^2MSE=n1∑i=1n∥x^pre,i−xpre,i∥22
(4) 结果分析根据多元回归模型预测的风化前化学成分，可以得到各文物样品在风化前的化学成分含量。通过与PCA分析的结果进行对比，验证预测结果的合理性。
python代码实现： import numpy as npimport pandas as pdfrom sklearn.decomposition import PCAfrom sklearn.linear_model import LinearRegressionfrom sklearn.metrics import mean_squared_errorfrom sklearn.model_selection import train_test_split# 假设 'data.csv' 包含玻璃文物的化学成分数据# 'type'列表示玻璃类型, 'weathered'列表示是否风化 (1表示风化, 0表示未风化)# 每行是一个样本，列为化学成分比例data = pd.read_csv('data.csv')# 1. 数据预处理# 我们需要对化学成分比例数据进行归一化处理def normalize_data(df):chemical_columns = df.columns.difference(['type', 'weathered'])df[chemical_columns] = df[chemical_columns].div(df[chemical_columns].sum(axis=1), axis=0)return dfdata_normalized = normalize_data(data)# 2. 主成分分析 (PCA)def perform_pca(df, n_components=3):chemical_columns = df.columns.difference(['type', 'weathered'])pca = PCA(n_components=n_components)pca_result = pca.fit_transform(df[chemical_columns])explained_variance = np.cumsum(pca.explained_variance_ratio_)print(f'Explained variance by PCA: {explained_variance}')# 返回降维后的数据和PCA对象return pd.DataFrame(pca_result, columns=[f'PC{i+1}' for i in range(n_components)]), pcapca_data, pca_model = perform_pca(data_normalized)# 将PCA结果添加回原数据data_normalized[['PC1', 'PC2', 'PC3']] = pca_data# 3. 多元线性回归模型# 我们将风化样本的成分作为输入，预测风化前的成分# 需要创建风化前和风化后的数据集# 假设化学成分从 'comp1', 'comp2', ..., 'compN' 是化学成分的列名chemical_columns = data_normalized.columns.difference(['type', 'weathered', 'PC1', 'PC2', 'PC3'])# 分割风化点（post-weathering）与未风化点（pre-weathering）post_weathering_data = data_normalized[data_normalized['weathered'] == 1][chemical_columns]pre_weathering_data = data_normalized[data_normalized['weathered'] == 0][chemical_columns]# 创建训练和测试集 (使用风化后的数据预测风化前的数据)X_train, X_test, y_train, y_test = train_test_split(post_weathering_data, pre_weathering_data, test_size=0.2, random_state=42)# 建立回归模型regressor = LinearRegression()regressor.fit(X_train, y_train)# 进行预测y_pred = regressor.predict(X_test)# 评估模型性能mse = mean_squared_error(y_test, y_pred)print(f'Mean Squared Error: {mse}')# 输出预测结果和实际结果对比comparison = pd.DataFrame({'Actual': y_test.sum(axis=1), 'Predicted': y_pred.sum(axis=1)})print(comparison.head())# 4. 结果分析与可视化import matplotlib.pyplot as plt# 可视化实际和预测结果对比plt.scatter(y_test.sum(axis=1), y_pred.sum(axis=1), alpha=0.5)plt.xlabel('Actual Pre-Weathering Composition')plt.ylabel('Predicted Pre-Weathering Composition')plt.title('Actual vs Predicted Pre-Weathering Composition')plt.show()问题2要求根据附件数据分析高钾玻璃和铅钡玻璃的分类规律，并对每个类别进行进一步的亚类划分。以下是详细的建模和解答过程，结合数学公式和复杂的分析方法来解决问题。
问题 2 解答过程分类规律分析： 我们首先要对高钾玻璃和铅钡玻璃进行初步分类分析，确定每类玻璃的化学成分特征。这可以通过对不同化学成分的分布进行统计描述，使用均值、方差等指标来识别两种玻璃类型之间的差异。
亚类划分： 对每类玻璃进一步进行亚类划分，可以借助聚类分析（如K-means、层次聚类）来识别在同一玻璃类型内部是否存在明显的分组。此外，我们可以使用判别分析（如线性判别分析，LDA）来确定哪些化学成分对类别划分最为重要。
敏感性分析： 最后，需要对分类和亚类划分结果的合理性和敏感性进行分析。通过调整不同的模型参数（如聚类中心个数、判别分析的阈值等），分析结果是否具有稳定性。
首先，对于原始数据，需要进行适当的处理。假设我们有化学成分的矩阵X∈Rn × pX \in \mathbb{R}^{n \times p}X∈Rn×p，其中nnn 是样本数，ppp 是化学成分的数量。每个样本的化学成分比例可能不完全精确，因此我们要进行归一化处理。
化学成分归一化公式为：
X i ′ = X i∑j=1pXij对于i=1,2,...,nX_i' = \frac{X_i}{\sum_{j=1}^{p} X_{ij}} \quad \text{对于} \quad i = 1, 2, ..., nXi′=∑j=1pXijXi对于i=1,2,...,n
其中，X i ′ X_i'Xi′ 表示第iii 个样本的归一化化学成分向量。
为了找出高钾玻璃和铅钡玻璃的分类规律，假设我们有两个类别C 1 C_1C1（高钾玻璃）和C 2 C_2C2（铅钡玻璃），对于每个化学成分X j X_jXj 我们可以分别计算各类别的均值μ C1 , j\mu_{C_1, j}μC1,j 和μ C2 , j\mu_{C_2, j}μC2,j，以及方差σ C1 , j2 \sigma_{C_1, j}^2σC1,j2 和σ C2 , j2 \sigma_{C_2, j}^2σC2,j2。
均值差异公式：
Δμ j =μ C1 , j−μ C2 , j\Delta \mu_j = \mu_{C_1, j} - \mu_{C_2, j}Δμj=μC1,j−μC2,j
我们可以通过计算每个化学成分的均值差异Δμ j \Delta \mu_jΔμj 来找出最具区分度的化学成分。如果∣Δμ j ∣|\Delta \mu_j|∣Δμj∣ 较大，说明该化学成分对于分类具有显著的作用。
为了进一步对高钾玻璃和铅钡玻璃进行亚类划分，我们可以使用K-means聚类，其目标是将同类玻璃样本进一步划分成kkk 个亚类。K-means的目标函数为：
J=∑i = 1k ∑x ∈Ci∥x−μ i ∥ 2 J = \sum_{i=1}^{k} \sum_{x \in C_i} \| x - \mu_i \|^2J=∑i=1k∑x∈Ci∥x−μi∥2
其中，μ i \mu_iμi 是第iii 个簇的中心，x∈C i x \in C_ix∈Ci 表示属于第iii 个簇的样本。我们通过迭代更新簇中心和样本分配，最小化目标函数JJJ。
对于每一类玻璃，我们可以分别选择不同的聚类数kkk，例如，高钾玻璃可能有两种亚类（如亚高钾和高钾），而铅钡玻璃可能有三种亚类。
为了找到哪些化学成分对类别划分最重要，可以使用线性判别分析（LDA）。LDA通过寻找一个投影向量www，使得不同类别之间的类间方差和类内方差的比值最大化：
LDA目标函数：
w=arg max w wTSB wwTSW ww = \arg \max_w \frac{w^T S_B w}{w^T S_W w}w=argmaxwwTSWwwTSBw
其中，S B S_BSB 是类间散布矩阵，S W S_WSW 是类内散布矩阵。它们的定义如下：
S B =∑i = 1k n i (μ i −μ)(μ i −μ) T S_B = \sum_{i=1}^{k} n_i (\mu_i - \mu)(\mu_i - \mu)^TSB=∑i=1kni(μi−μ)(μi−μ)T
S W =∑i = 1k ∑x ∈Ci(x−μ i )(x−μ i ) T S_W = \sum_{i=1}^{k} \sum_{x \in C_i} (x - \mu_i)(x - \mu_i)^TSW=∑i=1k∑x∈Ci(x−μi)(x−μi)T
通过最大化类间方差与类内方差的比值，LDA能找到化学成分中对分类最有效的方向。
六、分类结果的敏感性分析为了验证分类结果的合理性，我们可以进行交叉验证或留一法验证，从而评估模型在不同数据划分下的稳定性。定义分类准确率AAA 如下：
A=正确分类的样本数 总样本数 A = \frac{\text{正确分类的样本数}}{\text{总样本数}}A=总样本数正确分类的样本数
同时，通过调整聚类数kkk 或判别分析中的投影维度ddd，分析模型对参数变化的敏感性。
python代码实现： import pandas as pdfrom sklearn.preprocessing import StandardScaler# 假设我们已经有数据文件表单1和表单2# 数据格式示例 (csv 格式):# 表单 1: 包含文物的基本信息# 表单 2: 包含化学成分的比例信息df_glass_info = pd.read_csv('form1.csv')df_composition = pd.read_csv('form2.csv')# 填充缺失值为0，代表未检测到该成分df_composition.fillna(0, inplace=True)# 归一化每一行的化学成分，使总和为100df_composition_normalized = df_composition.div(df_composition.sum(axis=1), axis=0) * 100# 仅提取有用的化学成分列X = df_composition_normalized.values# 标准化化学成分特征scaler = StandardScaler()X_scaled = scaler.fit_transform(X)# 根据表单1中的类别信息进行初步分类# 假设 'category' 列代表玻璃类型，值为 '高钾玻璃' 或 '铅钡玻璃'df_glass_info['category'] = df_glass_info['category'].map({'高钾玻璃': 0, '铅钡玻璃': 1})# 分别提取高钾玻璃和铅钡玻璃的数据high_potassium_glass = df_composition_normalized[df_glass_info['category'] == 0]lead_barium_glass = df_composition_normalized[df_glass_info['category'] == 1]# 计算每种玻璃类型的均值和标准差mean_high_potassium = high_potassium_glass.mean()mean_lead_barium = lead_barium_glass.mean()std_high_potassium = high_potassium_glass.std()std_lead_barium = lead_barium_glass.std()# 输出化学成分的均值和标准差差异diff_means = mean_high_potassium - mean_lead_bariumprint("化学成分均值差异:\n", diff_means)from sklearn.cluster import KMeansimport numpy as np# 对高钾玻璃和铅钡玻璃分别进行聚类分析# 高钾玻璃分为2个亚类，铅钡玻璃分为3个亚类kmeans_high_potassium = KMeans(n_clusters=2, random_state=0)kmeans_lead_barium = KMeans(n_clusters=3, random_state=0)# 聚类结果high_potassium_labels = kmeans_high_potassium.fit_predict(high_potassium_glass)lead_barium_labels = kmeans_lead_barium.fit_predict(lead_barium_glass)print("高钾玻璃的聚类标签:\n", high_potassium_labels)print("铅钡玻璃的聚类标签:\n", lead_barium_labels)from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA# LDA 需要类别标签和化学成分数据X = df_composition_normalized.valuesy = df_glass_info['category'].values# 线性判别分析模型lda = LDA(n_components=1)X_lda = lda.fit_transform(X, y)print("LDA 投影后的数据:\n", X_lda)print("判别分析的线性系数:\n", lda.coef_)from sklearn.model_selection import cross_val_scorefrom sklearn.ensemble import RandomForestClassifier# 使用随机森林模型进行交叉验证clf = RandomForestClassifier(random_state=0)scores = cross_val_score(clf, X_scaled, y, cv=5)# 输出分类准确率print(f"交叉验证分类准确率: {np.mean(scores):.2f}")查看完整思路详见： 【腾讯文档】2024高教社杯全国大学生数学建模竞赛全题目深度解析（建模过程+代码实现+论文指导） https://docs.qq.com/doc/DSGdreXpIYlN2RUlZ