PCA_df = pd.DataFrame(PCAPipeline.fit_transform(X_train)) y_train.reset_index(drop=True, inplace=True) PCA_df = pd.concat([PCA_df, y_train], axis=1, ignore_index=True ) PCA_df.head(); plt.figure(figsize=(8,8)) sns.scatterplot(PCA_df[0],PCA_df[1],hue=PCA_df[3],palette=sns.color_palette("tab10", 2)) plt.show()
时间: 2024-03-18 17:44:13 浏览: 126
这段代码看起来像是使用PCA对训练数据进行降维,然后将降维后的数据和对应的标签合并成一个新的DataFrame,并画出散点图。其中,PCA_df是降维后的数据,X_train是原始数据,y_train是对应的标签,PCAPipeline是PCA模型的Pipeline。第一行代码使用fit_transform函数将训练数据进行降维,第二行代码将标签的索引重置并与降维后的数据进行合并,第三行代码画出散点图。
相关问题
import os import pandas as pd import joblib import numpy as np # 定义路径 input_dir = r'F:\Jungger 2D\test' # 新文件路径 model_dir = r'F:\Jungger 2D\data\PCA2\models' # 模型路径 output_dir = r'F:\Jungger 2D\data\PCA_results' # 输出路径 os.makedirs(output_dir, exist_ok=True) # 定义特征列 features = [ 'Heat Flow', 'erosion', 'Heat Flow232', 'Heat Flow66', 'Heat Flow145', 'Heat Flow298.9', 'Heat Flow23' ] # 加载模型 scaler = joblib.load(os.path.join(model_dir, 'global_scaler.pkl')) pca = joblib.load(os.path.join(model_dir, 'global_pca.pkl')) # 获取文件列表 files = [f for f in os.listdir(input_dir) if f.endswith('.txt')] # 处理每个文件 for file in files: file_path = os.path.join(input_dir, file) # 获取文件的完整路径 df = pd.read_csv(file_path, sep='\t') original_cols = df.iloc[:, :2] feature_data = df[features] # 显示前3行原始数据 print("原始特征数据表:") print(df[features].head(3)) scaled = scaler.transform(feature_data) # 打印标准化后的数据 print("标准化后的数据表:") print(scaled[:3]) # 打印前3行标准化后的数据 # 检查标准化后的数据的均值和标准差 print("标准化后的数据的均值:") print(np.mean(scaled, axis=0)) # 应接近0 print("标准化后的数据的标准差:") print(np.std(scaled, axis=0)) # 应接近1 # 执行PCA转换 pca_data = pca.transform(scaled) # 创建结果DataFrame pca_df = pd.DataFrame(pca_data, columns=[f'PC{i + 1}' for i in range(pca_data.shape[1])]) result_df = pd.concat([original_cols.reset_index(drop=True), pca_df.reset_index(drop=True)], axis=1) # 保存结果 output_file = os.path.join(output_dir, f"{file.split('.')[0]}_PCA.txt") try: result_df.to_csv(output_file, sep='\t', index=False) print(f"PCA结果已保存: {os.path.basename(output_file)}") except Exception as e: print(f"保存PCA结果文件 {output_file} 时发生错误:{e}") continue # 显示结果样例 print("结果样例:") print(result_df.head(3))这个是我用于pca的代码 标准化后的均值出现了100多
<think>好的,我现在需要解决用户在使用Python进行PCA时,经过StandardScaler标准化后的数据均值异常的问题,特别是当使用joblib加载scaler模型时,均值达到了100以上。这个问题看起来有点奇怪,因为标准化后的数据理论上均值应该接近0,标准差为1。我需要仔细分析可能的原因。
首先,用户提到他们使用了StandardScaler进行标准化,并且问题出现在加载保存的scaler模型时。可能的错误点可能在于数据的预处理步骤或者模型保存与加载的过程中出现了问题。
首先,StandardScaler的标准操作是减去均值并除以标准差,使得处理后的数据均值为0,标准差为1。所以如果标准化后的数据均值很大,可能说明标准化没有正确应用。可能的原因有几个:
1. **数据在标准化前未正确分割**:比如,用户可能在训练集和测试集划分之前进行了标准化,导致数据泄露,或者划分后未正确应用scaler。但用户提到是加载保存的scaler,因此可能是在预测时对新数据应用scaler时出了问题。
2. **保存和加载scaler时的版本或环境问题**:可能使用的joblib版本不兼容,导致scaler的参数没有被正确保存或加载。比如,不同版本的scikit-learn可能在保存模型时有差异,导致加载后的scaler参数错误。
3. **应用scaler时的数据形状不符**:例如,scaler是在二维数据上训练的(如形状为(n_samples, n_features)),但加载后应用的数据可能被错误地处理为一维数组,导致计算错误。
4. **手动修改了scaler的属性**:可能在保存scaler之前或之后,用户代码中意外修改了scaler的mean_或scale_属性,导致加载后的参数错误。
接下来,我需要验证这些可能性。例如,用户是否在训练时正确使用了fit_transform,而在测试时使用了transform?如果用户对测试数据错误地使用了fit_transform,会导致重新计算均值和标准差,而不是应用原来的参数,这可能引入错误。
另外,检查保存和加载scaler的代码是否正确。例如,保存时是否正确地使用joblib.dump,加载时是否用joblib.load。如果保存的scaler对象被破坏,或者加载时路径错误,可能导致scaler的均值和方差参数不正确。
另一个可能性是,用户在标准化后对数据进行了其他处理,例如错误的转换或计算,导致数据被改变。例如,是否在标准化之后又进行了某些操作,如乘以某个数值或添加了偏移量?
此外,检查数据的维度是否正确。例如,当使用scaler.transform时,输入的数据形状必须与训练时的数据形状一致。如果新数据的特征数量与训练时的不同,会导致错误,但此时通常会报错。不过,如果特征数量正确,但维度顺序错误(比如将行和列转置),也可能导致计算错误,从而使得标准化后的数据均值异常。
可能的解决步骤:
1. **检查scaler的参数是否正确保存和加载**:加载scaler后,打印scaler.mean_和scaler.scale_,确认其值与训练时的scaler是否一致。如果发现参数异常,可能是保存或加载过程中出了问题。
2. **验证标准化后的数据计算**:加载scaler后,手动对样本数据进行标准化计算,例如:(数据 - mean) / scale,确认结果是否与调用scaler.transform的结果一致。如果不一致,可能scaler的参数被篡改。
3. **检查数据预处理流程**:确保在训练和预测阶段都正确应用了transform方法,而不是fit或fit_transform。在新数据上错误调用fit_transform会导致scaler重新拟合数据,产生新的均值和标准差,这可能与原始数据分布不符,从而导致均值异常。
4. **确认数据形状和类型**:确保加载后的数据与训练时的数据具有相同的特征数和正确的形状。例如,输入的数据应该是二维数组,如(n_samples, n_features),避免一维数组导致的问题。
例如,假设用户保存了使用训练集X_train拟合的scaler,然后在加载后对新数据X_test应用transform。如果X_test的特征数或形状与X_train不同,就会导致错误。或者,在加载数据时,可能错误地将数据转换为错误的形状,比如转置了矩阵。
此外,检查是否有特征中存在大量缺失值或全为常数值的情况。如果某个特征在训练时的方差为0,StandardScaler会将其缩放为0,但在应用时如果该特征的方差不为0,可能会出现问题。但这种情况通常会导致警告或错误,而不是均值异常。
再比如,如果用户在保存scaler时使用了不同的库版本,例如scikit-learn版本升级导致参数保存格式变化,可能在加载时出现兼容性问题,导致scaler的参数错误。这时候需要确认训练环境和部署环境的库版本是否一致。
总结可能的原因:
- Scaler的参数在保存或加载过程中被破坏,导致均值和标准差错误。
- 在应用scaler时错误地调用了fit_transform而不是transform。
- 输入数据的形状或特征数与训练时不匹配。
- 数据预处理的其他步骤干扰了标准化结果。
- 环境或库版本不一致导致模型加载错误。
为了具体解决,用户应该逐步检查这些可能性。例如,在加载scaler后打印其属性,检查均值和标准差是否符合预期;检查新数据的处理流程是否正确,是否在正确的位置调用transform;确保数据形状正确,并且库版本一致。</think>### 问题分析:标准化后均值异常的可能原因
1. **错误的数据预处理流程**
- 若在加载`scaler`后对新数据使用`.fit_transform()`而非`.transform()`,会导致重新计算均值和标准差,而非应用原参数。此时若新数据分布与训练数据差异大,标准化后的均值可能异常[^3][^5]。
2. **模型保存与加载的兼容性问题**
- 使用`joblib`保存的`StandardScaler`模型,若加载时`scikit-learn`版本不一致,可能导致参数解析错误。例如,高版本保存的模型在低版本加载时,可能破坏`mean_`和`var_`属性[^4]。
3. **数据形状不匹配**
- 新数据输入维度与训练数据不一致(如特征数量不同),或未转换为二维数组(如形状为`(n_samples,)`而非`(n_samples, n_features)`),会导致计算错误[^2]。
4. **手动篡改`scaler`属性**
- 若代码中意外修改了`scaler.mean_`或`scaler.scale_`,加载后的参数将失效。
---
### 解决步骤与代码验证
#### 步骤1:检查`StandardScaler`参数
加载保存的`scaler`后,直接输出其属性:
```python
print("scaler.mean_:", scaler.mean_)
print("scaler.scale_:", scaler.scale_)
```
若输出值明显异常(如均值为100),说明模型保存/加载过程出错。
#### 步骤2:验证标准化计算逻辑
对样本数据手动计算标准化结果:
```python
import numpy as np
sample_data = np.array([[100, 200], [150, 250]]) # 示例数据
manual_result = (sample_data - scaler.mean_) / scaler.scale_
scaler_result = scaler.transform(sample_data)
print("手动计算结果:", manual_result)
print("scaler计算结果:", scaler_result)
```
若两者不一致,则`scaler`参数错误。
#### 步骤3:确保正确调用`transform`
检查代码中是否对新数据使用`.transform()`而非`.fit_transform()`:
```python
# 正确做法(加载后仅调用transform)
X_new_std = scaler.transform(X_new) # 不是fit_transform!
```
#### 步骤4:检查数据形状
确认输入数据为二维数组:
```python
# 错误示例:一维数据
X_wrong = [1, 2, 3]
# 正确示例:二维数据
X_correct = [[1, 2, 3], [4, 5, 6]]
```
#### 步骤5:版本一致性检查
确认训练与部署环境的库版本:
```bash
pip show scikit-learn # 检查版本是否一致
```
---
### 总结
若上述步骤均无问题,建议重新训练并保存`StandardScaler`模型,确保流程正确:
```python
from sklearn.preprocessing import StandardScaler
from joblib import dump, load
# 训练时保存
scaler = StandardScaler().fit(X_train)
dump(scaler, "scaler.joblib")
# 加载时验证
scaler_loaded = load("scaler.joblib")
X_test_std = scaler_loaded.transform(X_test)
print("标准化后均值:", X_test_std.mean(axis=0)) # 应接近0
```
---
# 设置中文支持 plt.rcParams['font.sans-serif'] = ['SimHei'] plt.rcParams['axes.unicode_minus'] = False mpl.rcParams['font.family'] = 'SimHei' 然后Excel有3个表单,表单1,2,3
<think>我们需要处理一个Excel文件,该文件包含3个工作表。我们将分别读取这三个工作表,然后进行数据预处理、分析和可视化。
由于数据具体内容未知,我们假设每个工作表的结构相似,包含相同类型的化学成分数据,可能还有分类信息(如类型、风化状态等)。
我们将按照以下步骤进行:
1. 读取Excel文件中的三个工作表,分别存储为三个DataFrame。
2. 对每个工作表进行数据预处理(处理缺失值、异常值等)。
3. 分析每个工作表的数据(描述性统计、成分分布、相关性等)。
4. 如果数据中包含分类信息(如玻璃类型、风化状态),则按类别进行分析。
5. 将三个工作表的数据合并,进行整体分析。
6. 可视化分析结果。
但是,由于问题描述中提到了“数据分类与归一化”,我们将重点放在数据预处理和特征工程上。
由于我们不知道具体列名,我们将假设每个工作表都有以下列(或类似):
- 类型(高钾、铅钡等)
- 风化状态(风化、未风化)
- 其他列为化学成分
如果实际列名不同,我们需要调整代码。
步骤:
1. 读取三个工作表
2. 数据预处理:
a. 检查并处理缺失值
b. 处理异常值(使用IQR方法)
c. 分类数据编码(如果需要)
3. 数据归一化:按类型和风化状态分组进行归一化(MinMax归一化或标准化)
4. 分别分析每个工作表,然后合并分析
由于数据量可能很大,我们将采用函数式编程,为重复操作编写函数。
注意:由于我们不知道具体列名,我们将通过特征(包含括号,如SiO2(%))来识别化学成分列。
我们将编写一个函数来处理每个工作表,然后进行合并分析。
但是,由于问题要求从数据分类与归一化开始,我们将重点放在这些步骤。
由于前面已经给出了一个完整的数据分析框架,这里我们针对三个工作表的情况进行调整。
由于时间关系,我们只给出整体框架,具体细节可以根据需要调整。
以下是完整代码:</think>下面是一个完整的解决方案,用于处理Excel文件中的三个表单数据。我将提供数据加载、预处理、分析和可视化的完整代码,并确保中文显示支持。
```python
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import matplotlib as mpl
from sklearn.preprocessing import MinMaxScaler, StandardScaler
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans
from scipy import stats
import re
import os
import warnings
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import classification_report, confusion_matrix
import joblib
# 设置中文支持
plt.rcParams['font.sans-serif'] = ['SimHei'] # 使用黑体
plt.rcParams['axes.unicode_minus'] = False # 解决负号显示问题
mpl.rcParams['font.family'] = 'SimHei' # 设置全局字体
sns.set(font='SimHei') # 设置Seaborn字体
# 忽略警告
warnings.filterwarnings('ignore')
def validate_column(df, column_name, possible_names):
"""验证列是否存在,支持模糊匹配"""
# 精确匹配
if column_name in df.columns:
return column_name
# 模糊匹配(忽略大小写、空格、下划线)
pattern = re.compile(r'[\s_]*'.join(possible_names), re.IGNORECASE)
for col in df.columns:
if pattern.search(col):
return col
# 列出所有可能的列名供参考
suggested = [col for col in df.columns if any(name.lower() in col.lower() for name in possible_names)]
raise ValueError(
f"缺失必要列:'{column_name}'\n"
f"可能的匹配列:{suggested}\n"
f"数据集所有列:{df.columns.tolist()}"
)
def load_all_sheets(file_path):
"""加载Excel文件中的所有工作表"""
print(f"加载Excel文件: {file_path}")
sheets = {}
try:
# 获取所有工作表名称
sheet_names = pd.ExcelFile(file_path).sheet_names
print(f"找到的工作表: {sheet_names}")
# 加载每个工作表
for sheet in sheet_names:
df = pd.read_excel(file_path, sheet_name=sheet)
sheets[sheet] = df
print(f"工作表 '{sheet}' 加载成功: {df.shape[0]}行, {df.shape[1]}列")
return sheets
except Exception as e:
raise ValueError(f"文件加载失败: {str(e)}")
def preprocess_sheet(df, sheet_name):
"""预处理单个工作表数据"""
print(f"\n预处理工作表: {sheet_name}")
# 识别关键列
try:
glass_type_col = validate_column(df, '类型', ['类型', '玻璃类型', 'glass', 'type', '玻璃种类'])
weathering_col = validate_column(df, '风化状态', ['风化', '风化状态', 'weathering', '状态', '风蚀情况'])
print(f"使用列名: 类型 -> '{glass_type_col}', 风化状态 -> '{weathering_col}'")
except ValueError as ve:
print(f"列验证错误: {str(ve)}")
glass_type_col, weathering_col = None, None
# 识别化学成分列
chem_cols = [col for col in df.columns if '(' in col and ')' in col]
if not chem_cols:
# 尝试其他方式识别化学成分列
chem_cols = [col for col in df.columns if any(keyword in col.lower() for keyword in ['成分', '含量', 'chem'])]
if not chem_cols:
chem_cols = [col for col in df.columns if col not in [glass_type_col, weathering_col] and df[col].dtype in ['float64', 'int64']]
print(f"警告: 无法准确识别化学成分列, 使用数值列代替: {chem_cols}")
print(f"识别到 {len(chem_cols)} 个化学成分列:")
print(chem_cols)
# 数据清洗
print("\n数据清洗...")
# 1. 处理缺失值 - 使用中位数填充
for col in chem_cols:
if df[col].isnull().sum() > 0:
median_val = df[col].median()
df[col].fillna(median_val, inplace=True)
print(f"列 '{col}' 使用中位数 {median_val:.4f} 填充了 {df[col].isnull().sum()} 个缺失值")
# 2. 处理异常值 - 使用IQR方法
for col in chem_cols:
q1 = df[col].quantile(0.25)
q3 = df[col].quantile(0.75)
iqr = q3 - q1
lower_bound = q1 - 1.5 * iqr
upper_bound = q3 + 1.5 * iqr
outliers = df[(df[col] < lower_bound) | (df[col] > upper_bound)]
if not outliers.empty:
# 使用中位数替换异常值
median_val = df[col].median()
df.loc[outliers.index, col] = median_val
print(f"列 '{col}' 修正了 {len(outliers)} 个异常值")
# 3. 确保风化状态列格式统一
if weathering_col and weathering_col in df.columns:
# 统一风化状态表示
df[weathering_col] = df[weathering_col].apply(
lambda x: '风化' if '风化' in str(x) else ('未风化' if '未风化' in str(x) else x)
)
print("\n风化状态分布:")
print(df[weathering_col].value_counts())
return df, chem_cols, glass_type_col, weathering_col
def classify_and_normalize(df, chem_cols, glass_type_col=None, weathering_col=None):
"""按类型和风化状态分类并进行归一化处理"""
print("\n执行分类与归一化处理...")
# 创建组合类型
if glass_type_col and weathering_col and glass_type_col in df.columns and weathering_col in df.columns:
df['组合类型'] = df[glass_type_col].astype(str) + '_' + df[weathering_col].astype(str)
elif glass_type_col and glass_type_col in df.columns:
df['组合类型'] = df[glass_type_col].astype(str)
elif weathering_col and weathering_col in df.columns:
df['组合类型'] = df[weathering_col].astype(str)
else:
df['组合类型'] = '所有样本'
print("警告: 未找到分类列,将所有样本视为同一类别")
# 归一化处理
scaler = MinMaxScaler()
normalized_cols = [f'归一化_{col}' for col in chem_cols]
# 按组合类型分组归一化
for category in df['组合类型'].unique():
mask = df['组合类型'] == category
temp = df.loc[mask, chem_cols]
# 检查是否有足够样本
if len(temp) < 2:
print(f"警告: 类别 '{category}' 样本不足 ({len(temp)}个), 跳过归一化")
continue
df.loc[mask, normalized_cols] = scaler.fit_transform(temp)
# 显示归一化后的统计信息
print("\n归一化后数据统计:")
display(df[normalized_cols].describe().T)
return df, chem_cols
def analyze_sheet(df, chem_cols, sheet_name, glass_type_col=None, weathering_col=None):
"""分析单个工作表数据"""
print(f"\n{'='*80}")
print(f"开始分析工作表: {sheet_name}")
# 1. 成分分布分析
plot_component_distributions(df, chem_cols, sheet_name)
# 2. 风化状态影响分析
if weathering_col and weathering_col in df.columns:
analyze_weathering_effects(df, chem_cols, weathering_col, sheet_name)
# 3. 玻璃类型分析
if glass_type_col and glass_type_col in df.columns:
analyze_glass_types(df, chem_cols, glass_type_col, weathering_col, sheet_name)
# 4. 相关性分析
analyze_correlations(df, chem_cols, sheet_name)
# 5. 聚类分析
df = perform_clustering(df, chem_cols, sheet_name)
# 6. PCA降维分析
perform_pca_analysis(df, chem_cols, sheet_name)
print(f"工作表 {sheet_name} 分析完成!\n")
return df
def plot_component_distributions(df, chem_cols, sheet_name):
"""绘制化学成分分布图"""
print(f"\n化学成分分布分析 - {sheet_name}:")
# 设置子图布局
n_cols = 4
n_rows = (len(chem_cols) + n_cols - 1) // n_cols
fig, axes = plt.subplots(n_rows, n_cols, figsize=(18, 4*n_rows))
axes = axes.flatten()
for i, col in enumerate(chem_cols):
ax = axes[i]
sns.histplot(df[col], kde=True, ax=ax, color='skyblue')
ax.set_title(f'{col}分布', fontsize=12)
ax.set_xlabel('')
# 隐藏多余的子图
for j in range(i+1, len(axes)):
axes[j].set_visible(False)
plt.suptitle(f'化学成分分布 - {sheet_name}', fontsize=16)
plt.tight_layout(rect=[0, 0, 1, 0.96]) # 为总标题留出空间
plt.savefig(f'{sheet_name}_化学成分分布.png', dpi=300)
plt.show()
# 箱线图汇总
plt.figure(figsize=(14, 8))
sns.boxplot(data=df[chem_cols], orient='h', palette='viridis')
plt.title(f'{sheet_name} - 化学成分分布箱线图', fontsize=16)
plt.tight_layout()
plt.savefig(f'{sheet_name}_化学成分箱线图.png', dpi=300)
plt.show()
def analyze_weathering_effects(df, chem_cols, weathering_col, sheet_name):
"""分析风化对化学成分的影响"""
print(f"\n风化状态对化学成分的影响分析 - {sheet_name}:")
# 按风化状态分组统计
weathering_stats = df.groupby(weathering_col)[chem_cols].agg(['mean', 'std'])
display(weathering_stats)
# 绘制风化状态对比
plt.figure(figsize=(14, 8))
melted_df = df.melt(id_vars=[weathering_col], value_vars=chem_cols,
var_name='成分', value_name='含量')
sns.boxplot(x='成分', y='含量', hue=weathering_col, data=melted_df, palette='Set2')
plt.title(f'{sheet_name} - 不同风化状态下化学成分含量对比', fontsize=16)
plt.xticks(rotation=45, ha='right')
plt.legend(title='风化状态')
plt.tight_layout()
plt.savefig(f'{sheet_name}_风化状态对比.png', dpi=300)
plt.show()
# 风化状态差异显著性检验
print("\n风化状态差异显著性检验 (独立样本t检验):")
results = []
for col in chem_cols:
unweathered = df[df[weathering_col] == '未风化'][col]
weathered = df[df[weathering_col] == '风化'][col]
if len(unweathered) < 2 or len(weathered) < 2:
print(f"警告: 列 '{col}' 样本不足,跳过检验")
continue
t_stat, p_value = stats.ttest_ind(unweathered, weathered, nan_policy='omit')
mean_diff = unweathered.mean() - weathered.mean()
results.append({
'成分': col,
'均值差(未风化-风化)': mean_diff,
't统计量': t_stat,
'p值': p_value,
'显著(p<0.05)': p_value < 0.05
})
display(pd.DataFrame(results))
def analyze_glass_types(df, chem_cols, glass_type_col, weathering_col, sheet_name):
"""分析不同玻璃类型的化学成分特征"""
print(f"\n玻璃类型化学成分特征分析 - {sheet_name}:")
# 按玻璃类型分组统计
type_stats = df.groupby(glass_type_col)[chem_cols].agg(['mean', 'std'])
display(type_stats)
# 可视化不同玻璃类型的化学成分
plt.figure(figsize=(16, 10))
melted_df = df.melt(id_vars=[glass_type_col], value_vars=chem_cols,
var_name='成分', value_name='含量')
sns.boxplot(x='成分', y='含量', hue=glass_type_col, data=melted_df, palette='tab10')
plt.title(f'{sheet_name} - 不同玻璃类型化学成分含量对比', fontsize=16)
plt.xticks(rotation=45, ha='right')
plt.legend(title='玻璃类型')
plt.tight_layout()
plt.savefig(f'{sheet_name}_玻璃类型对比.png', dpi=300)
plt.show()
# 玻璃类型与风化状态的交互分析
if weathering_col and weathering_col in df.columns:
print("\n玻璃类型与风化状态的交互分析:")
plt.figure(figsize=(16, 12))
g = sns.FacetGrid(df, col=glass_type_col, hue=weathering_col,
palette='Set1', col_wrap=3, height=4)
g.map(sns.scatterplot, chem_cols[0], chem_cols[1], alpha=0.7)
g.add_legend(title='风化状态')
plt.suptitle(f'{sheet_name} - 玻璃类型与风化状态的化学成分关系', y=1.02)
plt.tight_layout()
plt.savefig(f'{sheet_name}_玻璃类型与风化状态.png', dpi=300)
plt.show()
def analyze_correlations(df, chem_cols, sheet_name):
"""分析化学成分之间的相关性"""
print(f"\n化学成分相关性分析 - {sheet_name}:")
# 计算相关系数矩阵
corr_matrix = df[chem_cols].corr()
# 绘制热力图
plt.figure(figsize=(14, 12))
mask = np.triu(np.ones_like(corr_matrix, dtype=bool))
sns.heatmap(corr_matrix, mask=mask, annot=True, fmt=".2f",
cmap='coolwarm', center=0, vmin=-1, vmax=1)
plt.title(f'{sheet_name} - 化学成分相关性热力图', fontsize=16)
plt.tight_layout()
plt.savefig(f'{sheet_name}_化学成分相关性.png', dpi=300)
plt.show()
# 找出强相关关系
strong_correlations = corr_matrix.unstack().sort_values(ascending=False)
strong_correlations = strong_correlations[strong_correlations != 1] # 排除自身相关
strong_correlations = strong_correlations[abs(strong_correlations) > 0.7] # 只显示强相关
print("\n强相关关系 (|r| > 0.7):")
display(pd.DataFrame(strong_correlations).reset_index().rename(
columns={'level_0': '成分1', 'level_1': '成分2', 0: '相关系数'}))
def perform_clustering(df, chem_cols, sheet_name):
"""使用K-means进行聚类分析"""
print(f"\n聚类分析 (K-means) - {sheet_name}:")
# 数据标准化
scaler = StandardScaler()
scaled_data = scaler.fit_transform(df[chem_cols])
# 使用肘部法则确定最佳聚类数
distortions = []
K_range = range(2, 11)
for k in K_range:
kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)
kmeans.fit(scaled_data)
distortions.append(kmeans.inertia_)
plt.figure(figsize=(10, 6))
plt.plot(K_range, distortions, 'bo-')
plt.xlabel('聚类数 (k)')
plt.ylabel('距离平方和')
plt.title(f'{sheet_name} - 肘部法则确定最佳聚类数')
plt.grid(True)
plt.tight_layout()
plt.savefig(f'{sheet_name}_肘部法则.png', dpi=300)
plt.show()
# 选择聚类数(这里取肘部拐点,实际分析中应结合业务)
optimal_k = 4
kmeans = KMeans(n_clusters=optimal_k, random_state=42, n_init=10)
clusters = kmeans.fit_predict(scaled_data)
df['Cluster'] = clusters
# 聚类结果可视化
plt.figure(figsize=(10, 6))
if len(chem_cols) >= 2:
sns.scatterplot(x=chem_cols[0], y=chem_cols[1], hue='Cluster',
data=df, palette='viridis', s=100)
plt.title(f'{sheet_name} - K-means聚类结果 (k={optimal_k})', fontsize=16)
else:
# 如果只有一个成分,使用直方图
for cluster in range(optimal_k):
sns.histplot(df[df['Cluster'] == cluster][chem_cols[0]],
kde=True, label=f'Cluster {cluster}', alpha=0.5)
plt.title(f'{sheet_name} - K-means聚类结果 (k={optimal_k})', fontsize=16)
plt.legend()
plt.tight_layout()
plt.savefig(f'{sheet_name}_聚类结果.png', dpi=300)
plt.show()
# 聚类特征分析
print(f"\n各聚类化学成分特征 (k={optimal_k}):")
cluster_means = df.groupby('Cluster')[chem_cols].mean()
display(cluster_means)
return df
def perform_pca_analysis(df, chem_cols, sheet_name):
"""执行PCA降维分析"""
print(f"\n主成分分析 (PCA) - {sheet_name}:")
# 数据标准化
scaler = StandardScaler()
scaled_data = scaler.fit_transform(df[chem_cols])
# 执行PCA
pca = PCA()
pca_result = pca.fit_transform(scaled_data)
# 创建PCA结果DataFrame
pca_df = pd.DataFrame(data=pca_result, columns=[f'PC{i+1}' for i in range(pca_result.shape[1])])
pca_df = pd.concat([df.reset_index(drop=True), pca_df], axis=1)
# 解释方差比
plt.figure(figsize=(10, 6))
plt.bar(range(1, len(pca.explained_variance_ratio_)+1), pca.explained_variance_ratio_,
alpha=0.7, align='center', label='单个主成分方差')
plt.step(range(1, len(pca.explained_variance_ratio_)+1), np.cumsum(pca.explained_variance_ratio_),
where='mid', label='累计方差')
plt.ylabel('解释方差比')
plt.xlabel('主成分')
plt.title(f'{sheet_name} - 主成分解释方差')
plt.legend(loc='best')
plt.tight_layout()
plt.savefig(f'{sheet_name}_PCA方差解释率.png', dpi=300)
plt.show()
# 保留解释95%方差的主成分数
n_components = np.argmax(np.cumsum(pca.explained_variance_ratio_) >= 0.95) + 1
print(f"保留前 {n_components} 个主成分可解释 {np.cumsum(pca.explained_variance_ratio_)[n_components-1]:.2%} 的方差")
# 主成分载荷分析
loadings = pd.DataFrame(
pca.components_.T,
columns=[f'PC{i+1}' for i in range(len(chem_cols))],
index=chem_cols
)
print("\n前两个主成分的载荷:")
display(loadings[['PC1', 'PC2']])
# 绘制双标图
plt.figure(figsize=(12, 8))
if 'Cluster' in pca_df.columns:
sns.scatterplot(x='PC1', y='PC2', data=pca_df, hue='Cluster',
palette='viridis', s=100, alpha=0.7)
else:
sns.scatterplot(x='PC1', y='PC2', data=pca_df, s=100, alpha=0.7)
# 添加载荷向量
for i, feature in enumerate(chem_cols):
plt.arrow(0, 0, loadings.loc[feature, 'PC1']*3, loadings.loc[feature, 'PC2']*3,
color='r', alpha=0.5, head_width=0.05)
plt.text(loadings.loc[feature, 'PC1']*3.2, loadings.loc[feature, 'PC2']*3.2,
feature, color='darkred', fontsize=10)
plt.title(f'{sheet_name} - PCA双标图', fontsize=16)
plt.xlabel(f'PC1 ({pca.explained_variance_ratio_[0]:.2%})')
plt.ylabel(f'PC2 ({pca.explained_variance_ratio_[1]:.2%})')
plt.grid(True)
plt.tight_layout()
plt.savefig(f'{sheet_name}_PCA双标图.png', dpi=300)
plt.show()
def build_classification_model(df, chem_cols, glass_type_col):
"""构建玻璃类型分类模型"""
print("\n构建玻璃类型分类模型...")
# 准备数据
X = df[chem_cols]
y = df[glass_type_col]
# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.2, random_state=42, stratify=y
)
print(f"训练集大小: {X_train.shape[0]}, 测试集大小: {X_test.shape[0]}")
# 训练随机森林分类器
model = RandomForestClassifier(n_estimators=100, random_state=42)
model.fit(X_train, y_train)
# 评估模型
y_pred = model.predict(X_test)
print("\n分类报告:")
print(classification_report(y_test, y_pred))
# 混淆矩阵
plt.figure(figsize=(10, 8))
cm = confusion_matrix(y_test, y_pred)
sns.heatmap(cm, annot=True, fmt='d', cmap='Blues',
xticklabels=model.classes_, yticklabels=model.classes_)
plt.xlabel('预测值')
plt.ylabel('真实值')
plt.title('混淆矩阵')
plt.tight_layout()
plt.savefig('混淆矩阵.png', dpi=300)
plt.show()
# 特征重要性
feature_importances = pd.DataFrame({
'特征': chem_cols,
'重要性': model.feature_importances_
}).sort_values('重要性', ascending=False)
plt.figure(figsize=(12, 6))
sns.barplot(x='重要性', y='特征', data=feature_importances, palette='viridis')
plt.title('特征重要性排序')
plt.tight_layout()
plt.savefig('特征重要性.png', dpi=300)
plt.show()
# 保存模型
joblib.dump(model, 'glass_classifier.pkl')
print("模型已保存为 'glass_classifier.pkl'")
return model
# 主程序
if __name__ == "__main__":
# 原始数据路径
file_path = r"D:\Users\86157\Desktop\数学建模\附件.xlsx"
try:
# 1. 加载所有工作表
print("="*80)
print("【1/4】加载所有工作表...")
sheets = load_all_sheets(file_path)
processed_sheets = {}
# 2. 预处理每个工作表
print("\n" + "="*80)
print("【2/4】预处理每个工作表...")
for sheet_name, df in sheets.items():
print(f"\n处理工作表: {sheet_name}")
df, chem_cols, glass_type_col, weathering_col = preprocess_sheet(df.copy(), sheet_name)
df, chem_cols = classify_and_normalize(df, chem_cols, glass_type_col, weathering_col)
processed_sheets[sheet_name] = (df, chem_cols, glass_type_col, weathering_col)
# 3. 分析每个工作表
print("\n" + "="*80)
print("【3/4】分析每个工作表...")
analysis_results = {}
for sheet_name, (df, chem_cols, glass_type_col, weathering_col) in processed_sheets.items():
df = analyze_sheet(df, chem_cols, sheet_name, glass_type_col, weathering_col)
analysis_results[sheet_name] = df
# 4. 合并所有工作表进行综合分析
print("\n" + "="*80)
print("【4/4】合并所有工作表进行综合分析...")
# 合并所有数据
all_data = pd.concat([df for df, _, _, _ in processed_sheets.values()], ignore_index=True)
print(f"合并后总数据量: {all_data.shape[0]}行, {all_data.shape[1]}列")
# 识别公共化学成分列
common_chem_cols = set.intersection(*[set(chem_cols) for _, chem_cols, _, _ in processed_sheets.values()])
common_chem_cols = sorted(common_chem_cols)
print(f"共有 {len(common_chem_cols)} 个公共化学成分列:")
print(common_chem_cols)
# 分析合并后的数据
glass_type_col = validate_column(all_data, '类型', ['类型', '玻璃类型', 'glass', 'type', '玻璃种类'])
weathering_col = validate_column(all_data, '风化状态', ['风化', '风化状态', 'weathering', '状态', '风蚀情况'])
all_data, _ = classify_and_normalize(all_data, common_chem_cols, glass_type_col, weathering_col)
all_data = analyze_sheet(all_data, common_chem_cols, "合并数据", glass_type_col, weathering_col)
# 构建分类模型
if glass_type_col and glass_type_col in all_data.columns:
build_classification_model(all_data, common_chem_cols, glass_type_col)
# 保存处理后的数据
print("\n保存处理后的数据...")
with pd.ExcelWriter('处理后的数据.xlsx') as writer:
for sheet_name, (df, _, _, _) in processed_sheets.items():
df.to_excel(writer, sheet_name=sheet_name, index=False)
all_data.to_excel(writer, sheet_name='合并数据', index=False)
print("处理后的数据已保存至 '处理后的数据.xlsx'")
print("\n所有处理和分析完成!")
except Exception as e:
print(f"\n{'*'*50}")
print(f"程序执行出错: {str(e)}")
print(f"{'*'*50}")
import traceback
traceback.print_exc()
```
### 代码功能详解
#### 1. 中文显示支持
```python
# 设置中文支持
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
mpl.rcParams['font.family'] = 'SimHei'
sns.set(font='SimHei')
```
这部分代码确保所有图表和标签都能正确显示中文,解决了matplotlib和seaborn中文显示问题。
#### 2. 多工作表处理
```python
def load_all_sheets(file_path):
"""加载Excel文件中的所有工作表"""
sheets = {}
sheet_names = pd.ExcelFile(file_path).sheet_names
for sheet in sheet_names:
df = pd.read_excel(file_path, sheet_name=sheet)
sheets[sheet] = df
return sheets
```
此函数加载Excel文件中的所有工作表,并返回一个字典,键为工作表名,值为对应的DataFrame。
#### 3. 智能列识别
```python
def validate_column(df, column_name, possible_names):
"""验证列是否存在,支持模糊匹配"""
# 模糊匹配逻辑...
```
该函数智能识别关键列(如"类型"、"风化状态"),支持多种常见列名变体,确保代码对不同格式数据的兼容性。
#### 4. 数据预处理流程
1. **缺失值处理**:使用中位数填充缺失值
2. **异常值处理**:基于IQR方法检测并修正异常值
3. **数据归一化**:按玻璃类型和风化状态分组进行MinMax归一化
4. **数据清洗**:统一风化状态的表示("风化"/"未风化")
#### 5. 全面分析模块
每个工作表独立分析以下内容:
- **成分分布**:直方图和箱线图展示
- **风化影响**:风化/未风化状态对比及t检验
- **玻璃类型**:不同类型玻璃的成分特征
- **相关性**:化学成分间相关性热力图
- **聚类分析**:K-means聚类及可视化
- **PCA降维**:主成分分析及双标图
#### 6. 多表合并分析
将所有工作表数据合并后进行:
- 公共化学成分分析
- 整体数据分布分析
- 玻璃类型分类模型构建(随机森林)
- 模型评估及特征重要性分析
### 输出结果
程序将为每个工作表和合并数据生成以下文件:
1. `{sheet_name}_化学成分分布.png` - 成分分布直方图
2. `{sheet_name}_化学成分箱线图.png` - 成分分布箱线图
3. `{sheet_name}_风化状态对比.png` - 风化状态对比图
4. `{sheet_name}_玻璃类型对比.png` - 玻璃类型对比图
5. `{sheet_name}_化学成分相关性.png` - 相关性热力图
6. `{sheet_name}_肘部法则.png` - K-means聚类肘部法则图
7. `{sheet_name}_聚类结果.png` - 聚类可视化
8. `{sheet_name}_PCA方差解释率.png` - PCA解释方差图
9. `{sheet_name}_PCA双标图.png` - PCA双标图
10. `处理后的数据.xlsx` - 预处理后的完整数据
11. `glass_classifier.pkl` - 玻璃类型分类模型
12. `混淆矩阵.png` - 模型混淆矩阵
13. `特征重要性.png` - 特征重要性排序
### 使用说明
1. 设置数据路径:
```python
file_path = r"D:\Users\86157\Desktop\数学建模\附件.xlsx"
```
2. 运行程序:
- 代码将自动处理所有工作表
- 每个工作表独立分析后合并分析
- 所有结果保存为图像和Excel文件
3. 自定义调整:
- 修改`classify_and_normalize`函数中的归一化方法
- 调整`perform_clustering`中的聚类数
- 修改`build_classification_model`中的分类器参数
此解决方案提供了从数据加载到最终分析的完整流程,特别针对多工作表Excel文件和玻璃文物数据的特点进行了优化。
阅读全文
相关推荐



