2024高教社杯全国大学生数学建模竞赛（C题）深度剖析 _ 建模完整过程+详细思路+代码全解析

问题 1 解答过程

建模方法：

为了解决该问题，采用主成分分析（PCA）和多元线性回归模型来进行建模。PCA用于降维，帮助我们识别与玻璃类型、风化状态、纹饰、颜色相关的主要成分，进而简化数据结构。多元线性回归则用于基于风化后的成分预测风化前的化学成分含量。

---

主成分分析（PCA）

(1) 数据预处理

首先，将每个样本的化学成分比例矩阵定义为 $X\in {\mathbb{R}}^{n\times p}$，其中 $n$ 是样本数量，$p$ 是化学成分的数量。由于数据中的成分累加和并不总是100%，对数据进行归一化处理。假设第 $i$ 个样本的化学成分比例为 ${\mathbf{x}}_i=[x_{i1},x_{i2},\dots ,x_{ip}]$，则归一化后的成分表示为：

$x_{ij}^{\prime }=\frac{x_{ij}}{{\sum }_{j=1}^p{x}_{ij}}$

确保成分比例的总和接近1。

(2) 协方差矩阵计算

在完成数据归一化后，构建数据矩阵 $X^{\prime }$ 并计算其协方差矩阵 $\Sigma$：

$\Sigma =\frac{1}{n-1}{\sum }_{i=1}^n({\mathbf{x}}_i^{\prime }-{\bar{\mathbf{x}}}^{\prime })({\mathbf{x}}_i^{\prime }-{\bar{\mathbf{x}}}^{\prime }{)}^{\top }$

其中，${\bar{\mathbf{x}}}^{\prime }$ 是样本的均值向量。协方差矩阵 $\Sigma$ 用于捕捉成分之间的相关性。

(3) 特征值分解

接下来，进行特征值分解（EVD），即：

$\Sigma =V\Lambda V^{\top }$

其中，$\Lambda =\text{diag}({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_p)$ 为协方差矩阵的特征值，${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_p$，$V=[{\mathbf{v}}_1,{\mathbf{v}}_2,\dots ,{\mathbf{v}}_p]$ 是协方差矩阵的特征向量。特征向量 ${\mathbf{v}}_i$ 对应主成分方向，而特征值 ${\lambda }_i$ 则表明该主成分解释的方差大小。

(4) 主成分选择

为了简化模型，只选择前 $k$ 个特征值较大的主成分。累积贡献率表示为：

$\text{累积贡献率}=\frac{{\sum }_{i=1}^k{\lambda }_i}{{\sum }_{i=1}^p{\lambda }_i}$

通常选择使累积贡献率达到80%-90%的 $k$ 值。选定的主成分矩阵为 $Z=X{V}_k$，其中 $V_k$ 是前 $k$ 个主成分方向组成的矩阵。

(5) PCA结果分析

将PCA降维后的数据 $Z$ 与玻璃类型、颜色、纹饰进行关联分析。使用PCA降维后的二维或三维散点图观察不同玻璃类型（如高钾玻璃、铅钡玻璃）的分布，以及不同风化程度与颜色的关联性。

---

多元线性回归模型

(1) 模型设定

我们要基于风化后的化学成分预测其风化前的成分。假设风化后的成分为 ${\mathbf{x}}_{\text{post}}=[x_{\text{post},1},x_{\text{post},2},\dots ,x_{\text{post},p}]$，风化前的成分为 ${\mathbf{x}}_{\text{pre}}=[x_{\text{pre},1},x_{\text{pre},2},\dots ,x_{\text{pre},p}]$，则可以构建多元线性回归模型：

${\mathbf{x}}_{\text{pre}}=\mathbf{B}{\mathbf{x}}_{\text{post}}+ϵ$

其中，$\mathbf{B}\in {\mathbb{R}}^{p\times p}$ 为回归系数矩阵，$ϵ$ 为残差项，表示模型误差。

(2) 回归系数的估计

为了估计回归系数矩阵 $\mathbf{B}$，采用最小二乘法（OLS, Ordinary Least Squares），即：

$\hat{\mathbf{B}}=\arg {\min }_{\mathbf{B}}{\sum }_{i=1}^n\Vert {\mathbf{x}}_{\text{pre},i}-\mathbf{B}{\mathbf{x}}_{\text{post},i}{\Vert }_2^2$

其解析解为：

$\hat{\mathbf{B}}=({\mathbf{X}}_{\text{post}}^{\top }{\mathbf{X}}_{\text{post}}{)}^{-1}{\mathbf{X}}_{\text{post}}^{\top }{\mathbf{X}}_{\text{pre}}$

其中，${\mathbf{X}}_{\text{post}}$ 和 ${\mathbf{X}}_{\text{pre}}$ 分别为风化后和风化前成分的数据矩阵。

(3) 模型评估与预测

通过交叉验证评估模型的预测能力。将模型应用于风化后的新数据，预测其风化前的化学成分：

${\hat{\mathbf{x}}}_{\text{pre}}=\hat{\mathbf{B}}{\mathbf{x}}_{\text{post}}$

对比实际数据与预测结果，计算预测误差，例如使用均方误差（MSE）评估：

$\text{MSE}=\frac{1}{n}{\sum }_{i=1}^n\Vert {\hat{\mathbf{x}}}_{\text{pre},i}-{\mathbf{x}}_{\text{pre},i}{\Vert }_2^2$

(4) 结果分析

根据多元回归模型预测的风化前化学成分，可以得到各文物样品在风化前的化学成分含量。通过与PCA分析的结果进行对比，验证预测结果的合理性。

python代码实现：

import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
from 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 df

data_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)]), pca

pca_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 解答过程

1. 分类规律分析：
   我们首先要对高钾玻璃和铅钡玻璃进行初步分类分析，确定每类玻璃的化学成分特征。这可以通过对不同化学成分的分布进行统计描述，使用均值、方差等指标来识别两种玻璃类型之间的差异。

2. 亚类划分：
   对每类玻璃进一步进行亚类划分，可以借助聚类分析（如K-means、层次聚类）来识别在同一玻璃类型内部是否存在明显的分组。此外，我们可以使用判别分析（如线性判别分析，LDA）来确定哪些化学成分对类别划分最为重要。

3. 敏感性分析：
   最后，需要对分类和亚类划分结果的合理性和敏感性进行分析。通过调整不同的模型参数（如聚类中心个数、判别分析的阈值等），分析结果是否具有稳定性。

首先，对于原始数据，需要进行适当的处理。假设我们有化学成分的矩阵 $X\in {\mathbb{R}}^{n\times p}$，其中 $n$ 是样本数，$p$ 是化学成分的数量。每个样本的化学成分比例可能不完全精确，因此我们要进行归一化处理。

化学成分归一化公式为：

$X_i^{\prime }=\frac{X_i}{{\sum }_{j=1}^p{X}_{ij}}\text{对于}i=1,2,...,n$

其中，$X_i^{\prime }$ 表示第 $i$ 个样本的归一化化学成分向量。

为了找出高钾玻璃和铅钡玻璃的分类规律，假设我们有两个类别 $C_1$（高钾玻璃）和 $C_2$（铅钡玻璃），对于每个化学成分 $X_j$ 我们可以分别计算各类别的均值 ${\mu }_{C_1,j}$ 和 ${\mu }_{C_2,j}$，以及方差 ${\sigma }_{C_1,j}^2$ 和 ${\sigma }_{C_2,j}^2$。

均值差异公式：

$\Delta {\mu }_j={\mu }_{C_1,j}-{\mu }_{C_2,j}$

我们可以通过计算每个化学成分的均值差异 $\Delta {\mu }_j$ 来找出最具区分度的化学成分。如果 $|\Delta {\mu }_j|$ 较大，说明该化学成分对于分类具有显著的作用。

为了进一步对高钾玻璃和铅钡玻璃进行亚类划分，我们可以使用K-means聚类，其目标是将同类玻璃样本进一步划分成 $k$ 个亚类。K-means的目标函数为：

$J={\sum }_{i=1}^k{\sum }_{x\in C_i}\Vert x-{\mu }_i{\Vert }^2$

其中，${\mu }_i$ 是第 $i$ 个簇的中心，$x\in C_i$ 表示属于第 $i$ 个簇的样本。我们通过迭代更新簇中心和样本分配，最小化目标函数 $J$。

对于每一类玻璃，我们可以分别选择不同的聚类数 $k$，例如，高钾玻璃可能有两种亚类（如亚高钾和高钾），而铅钡玻璃可能有三种亚类。

为了找到哪些化学成分对类别划分最重要，可以使用线性判别分析（LDA）。LDA通过寻找一个投影向量 $w$，使得不同类别之间的类间方差和类内方差的比值最大化：

LDA目标函数：

$w=\arg {\max }_w\frac{w^T{S}_{B}w}{w^T{S}_{W}w}$

其中，$S_B$ 是类间散布矩阵，$S_W$ 是类内散布矩阵。它们的定义如下：

$S_B={\sum }_{i=1}^k{n}_i({\mu }_i-\mu )({\mu }_i-\mu {)}^T$

$S_W={\sum }_{i=1}^k{\sum }_{x\in C_i}(x-{\mu }_i)(x-{\mu }_i{)}^T$

通过最大化类间方差与类内方差的比值，LDA能找到化学成分中对分类最有效的方向。

六、分类结果的敏感性分析

为了验证分类结果的合理性，我们可以进行交叉验证或留一法验证，从而评估模型在不同数据划分下的稳定性。定义分类准确率 $A$ 如下：

$A=\frac{\text{正确分类的样本数}}{\text{总样本数}}$

同时，通过调整聚类数 $k$ 或判别分析中的投影维度 $d$，分析模型对参数变化的敏感性。

python代码实现：

import pandas as pd
from 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)

# 归一化每一行的化学成分，使总和为100
df_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_barium
print("化学成分均值差异:\n", diff_means)

from sklearn.cluster import KMeans
import 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.values
y = 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_score
from 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}")

关注Unicorn建模，以往回答中有全部解答～