请给我2022年c题第二题的详细解题思路过程和代码

<think>我们正在讨论的是2022年全国大学生数学建模竞赛C题（玻璃文物风化分析）的第二问。根据用户需求，第二问包含两个部分： 1. 分析玻璃文物表面风化与其类型、纹饰和颜色的关系 2. 预测风化前的化学成分 根据之前提供的思路，我们将详细展开解题过程和代码实现。 注意：由于实际数据需要从官方提供的excel中读取，这里我们假设数据文件名为"glass_data.xlsx"，并包含以下字段： - 表面风化（有风化/无风化） - 玻璃类型（高钾/铅钡） - 纹饰（A类、B类等，具体看数据） - 颜色（黑色、浅蓝色等，具体看数据） - 化学成分（包括风化前后的各种化学成分含量，如二氧化硅(SiO2)、氧化钠(Na2O)等） 第二问的详细步骤： 第一部分：关系分析（风化与类型、纹饰、颜色） 步骤1：数据预处理 处理缺失值：对于分类变量（纹饰、颜色）用众数填充；对于数值变量（化学成分）用中位数或均值填充（根据分布决定）。 将分类变量进行编码（如纹饰和颜色可能有多个类别，使用独热编码或标签编码）。 步骤2：统计检验 a) 对于玻璃类型与风化的关系：因为玻璃类型只有两类（高钾、铅钡），使用卡方检验判断是否独立。 b) 对于纹饰与风化的关系：纹饰可能有多个类别，同样使用卡方检验，但要注意期望频数不能太小（如果太小，考虑合并类别或使用Fisher精确检验）。 c) 对于颜色与风化的关系：同上。 步骤3：计算关联强度 使用Cramer's V系数来量化关联程度（适用于大于2×2的列联表）。 步骤4：可视化 绘制柱状图或堆叠柱状图展示不同类别下的风化比例。 第二部分：化学成分预测（由风化后的成分预测风化前的成分） 注意：我们有一些样本同时有风化前后的化学成分数据，可以用于建模。 步骤1：数据准备 提取同时有风化和未风化成分的样本（即配对样本）。 将数据分为两组：高钾玻璃和铅钡玻璃（因为风化机制可能不同，建议分别建模）。 步骤2：建立预测模型 模型1：流失率模型（简单线性模型） 对于每个化学成分，计算未风化样本和配对的风化样本的差异，得到平均流失率（或流失比例）。 然后对于只有风化成分的样本，用公式：未风化成分 = 风化成分 / (1 - 流失率) 注意：流失率不能为负，且如果某种元素在风化过程中增加（如土壤中的元素渗入），则流失率可能为负。 模型2：多元线性回归（或岭回归） 输入：风化后的各化学成分含量（可能还有类型、纹饰、颜色等作为辅助变量） 输出：未风化的各化学成分含量 由于化学成分之间可能存在多重共线性，使用岭回归（Ridge Regression）可以避免过拟合。 步骤3：模型评估 使用配对样本的一部分作为训练集，另一部分作为测试集，计算均方误差（MSE）或平均绝对误差（MAE）。 步骤4：预测 对未知的样本（只有风化成分）进行预测。 代码实现： 第一部分：关系分析</think>### 2022年全国大学生数学建模竞赛C题第二问详细解析 **问题要求**：分析玻璃文物表面风化与类型、纹饰、颜色的关系，并预测风化前化学成分。以下是完整解题思路与代码实现： --- ### 一、关系分析（风化与类型/纹饰/颜色） #### 核心方法：统计检验 + 关联强度量化 1. **数据预处理** - 处理缺失值：分类变量用众数填充，数值变量用中位数填充 - 分类变量编码：纹饰/颜色转换为哑变量 ```python import pandas as pd import numpy as np from scipy.stats import chi2_contingency # 数据加载与预处理 data = pd.read_excel('glass_data.xlsx') data['颜色'] = data['颜色'].fillna(data['颜色'].mode()[0]) data['纹饰'] = data['纹饰'].fillna(data['纹饰'].mode()[0]) data = pd.get_dummies(data, columns=['纹饰', '颜色']) # 哑变量编码 ``` 2. **卡方检验与Cramer's V系数** $$ \chi^2 = \sum\frac{(O_{ij} - E_{ij})^2}{E_{ij}}, \quad V = \sqrt{\frac{\chi^2}{n \times \min(r-1,c-1)}} $$ - $O_{ij}$: 观测频数 - $E_{ij}$: 期望频数 - $n$: 样本量 - $r,c$: 行列维度数 ```python def cramers_v(conf_matrix): chi2 = chi2_contingency(conf_matrix)[0] n = conf_matrix.sum().sum() phi2 = chi2 / n r, k = conf_matrix.shape return np.sqrt(phi2 / min((k-1), (r-1))) # 分析玻璃类型与风化的关系 type_table = pd.crosstab(data['表面风化'], data['玻璃类型']) chi2, p, _, _ = chi2_contingency(type_table) v_type = cramers_v(type_table.values) print(f"类型关联强度: V={v_type:.3f} (p={p:.4f})") # 分析所有纹饰类别 for pattern in data.filter(like='纹饰_').columns: pattern_table = pd.crosstab(data['表面风化'], data[pattern]) if pattern_table.sum().min() > 5: # 确保最小期望频数>5 v_pattern = cramers_v(pattern_table.values) print(f"纹饰[{pattern}]关联强度: V={v_pattern:.3f}") ``` 3. **结果解读** - $V < 0.1$: 无关联 - $0.1 \leq V < 0.3$: 弱关联 - $0.3 \leq V < 0.5$: 中等关联 - $V \geq 0.5$: 强关联 - 示例输出： ``` 类型关联强度: V=0.372 (p=0.0002) → 中等关联 纹饰_几何纹关联强度: V=0.412 → 中等关联 ``` --- ### 二、化学成分预测（风化后→风化前） #### 核心方法：流失率模型 + 岭回归 1. **数据准备** - 提取配对样本（同时有风化和未风化数据） - 按玻璃类型分组建模（高钾/铅钡风化机制不同） ```python from sklearn.linear_model import Ridge from sklearn.preprocessing import StandardScaler # 提取配对样本 paired_data = data[data['风化状态']=='同时有'] X_weathered = paired_data.filter(regex='风化点_.*') # 风化成分 y_original = paired_data.filter(regex='未风化_.*') # 未风化成分 # 按玻璃类型分组 high_k = paired_data[paired_data['玻璃类型']=='高钾'] lead_barium = paired_data[paired_data['玻璃类型']=='铅钡'] ``` 2. **流失率模型** $$ C_{\text{未风化}} = \frac{C_{\text{风化点}}}{1 - \beta_i} $$ $\beta_i$为元素$i$的平均流失率： ```python # 计算各元素流失率 loss_rates = 1 - (X_weathered.mean() / y_original.mean()) print("元素平均流失率:\n", loss_rates) # 应用流失率预测 def predict_by_loss_rate(weathered_data): return weathered_data.div(1 - loss_rates, axis=1) ``` 3. **岭回归模型** $$ \mathbf{Y} = \mathbf{X}\beta + \epsilon $$ - $\mathbf{Y}$: 未风化成分矩阵 - $\mathbf{X}$: 风化成分矩阵 ```python # 数据标准化 scaler = StandardScaler() X_scaled = scaler.fit_transform(X_weathered) y_scaled = scaler.fit_transform(y_original) # 建立岭回归模型（L2正则化防止过拟合） model = Ridge(alpha=0.5) # 正则化强度 model.fit(X_scaled, y_scaled) # 预测新样本 new_sample = scaler.transform([[70.1, 12.3, ...]]) # 输入风化成分 predicted = model.predict(new_sample) original_comp = scaler.inverse_transform(predicted) print("预测未风化成分:", original_comp) ``` 4. **模型评估** ```python from sklearn.metrics import mean_absolute_error, r2_score from sklearn.model_selection import train_test_split # 划分训练集/测试集 X_train, X_test, y_train, y_test = train_test_split(X_scaled, y_scaled, test_size=0.2) model.fit(X_train, y_train) y_pred = model.predict(X_test) print("R² Score:", r2_score(y_test, y_pred)) print("MAE:", mean_absolute_error(y_test, y_pred)) ``` --- ### 三、关键解题要点 1. **数据特殊性处理** - 铅钡玻璃需单独建模：其$\ce{PbO}$含量在风化中显著变化 - 高钾玻璃关注$\ce{K2O}$、$\ce{SiO2}$流失率 2. **可视化验证** - 三元图分析$\ce{SiO2}-\ce{Na2O}-\ce{CaO}$风化轨迹 ```python import plotly.express as px fig = px.scatter_ternary(data, a="SiO2", b="Na2O", c="CaO", color="表面风化", symbol="玻璃类型") fig.show() ``` 3. **模型选择依据** | 模型 | 适用场景 | 优势 | |------------|-------------------------|----------------------| | 流失率模型 | 元素流失机制明确 | 物理意义清晰 | | 岭回归 | 多元素协同变化 | 捕捉元素间交互作用 | --- ### 四、完整代码结构 ```python # 1. 数据预处理 data = load_and_preprocess('glass_data.xlsx') # 2. 关系分析 analyze_relationship(data, ['玻璃类型', '纹饰_A', '纹饰_B', '颜色_绿']) # 3. 化学成分预测 # 高钾玻璃模型 high_k_model = train_ridge_model(high_k_data) # 铅钡玻璃模型 pb_model = train_ridge_model(lead_barium_data) # 4. 预测未知样本 predict_original_composition(new_sample, glass_type='铅钡') ``` > **注**：实际应用需根据官方数据调整特征选择，重点关注$\ce{SiO2}$、$\ce{PbO}$、$\ce{K2O}$等核心成分[^1][^2]。 --- ### 五、结果应用场景 1. 文物修复：预测原始成分指导复原材料配比 2. 文物鉴定：通过风化特征辅助断代 3. 保护策略：针对不同玻璃类型制定差异化保护方案