PCA算法降维代码示例

这段代码将数据进行PCA降维至3维，并绘制一个三维散点图，展示降维后的前3个主成分。

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
import matplotlib.colors as mcolors
from mpl_toolkits.mplot3d import Axes3D

# 读取数据
file_path = '4_SmCrTe3_Study_AFM_Select.txt'
data = pd.read_csv(file_path, sep='\t', header=None)

# 命名列
columns = ['ID', 'f1', 'f2', 'f3', 'f4', 'f5', 'f6', 'f7', 'Energy', 'Unused']
data.columns = columns

# 删除不需要的列
data = data.drop(columns=['ID', 'Unused'])

# 数据概览
print(data.describe())

# 分析Energy列的统计数据
energy_stats = data['Energy'].describe()
print("\nEnergy column statistics:")
print(energy_stats)

# 1. 直方图（1_Energy_Analysis_Histogram.png）
plt.figure(figsize=(12, 6))

# 直方图
plt.subplot(1, 2, 1)
sns.histplot(data['Energy'], kde=True)
plt.title('Energy Distribution')
plt.xlabel('Energy')

# 在直方图中标注count数量
for patch in plt.gca().patches:
    height = patch.get_height()
    plt.annotate(f'{height:.0f}', (patch.get_x() + patch.get_width() / 2, height), ha='center', va='bottom')

# 第二个直方图，用于替代箱线图
plt.subplot(1, 2, 2)
sns.histplot(data['Energy'], bins=30, kde=True)
plt.title('Energy Distribution (Detailed)')
plt.xlabel('Energy')

# 在直方图中标注count数量
for patch in plt.gca().patches:
    height = patch.get_height()
    plt.annotate(f'{height:.0f}', (patch.get_x() + patch.get_width() / 2, height), ha='center', va='bottom')

plt.tight_layout()
plt.show()

# 检查并处理NaN值
print("\nNumber of NaN values in each column:")
print(data.isna().sum())

# 使用插值方法填补NaN值
data = data.interpolate()

# 再次检查NaN值是否已经处理
print("\nNumber of NaN values in each column after interpolation:")
print(data.isna().sum())

# 2. 散点图（2_Energy_Analysis_Scatter.png）
plt.figure(figsize=(12, 6))
sns.scatterplot(data=data, x=data.index, y='Energy', color='dodgerblue')
plt.title('Selected SmCrTe3 Energy Distribution', fontsize=15)
plt.xlabel('Sample Index', fontsize=12)
plt.ylabel('Energy (meV)', fontsize=12)
plt.show()

# 3. 热力图（3_Single_f-Orbital_Couplings_with_Energy_Hot.png）
plt.figure(figsize=(12, 8))
sns.heatmap(data.corr(), annot=True, cmap='coolwarm', center=0, linewidths=0.5)
plt.title('Correlation Matrix of f-Orbital Occupations and Energy', fontsize=15)
plt.show()

# 双轨道和能量关系（4_Double_f-Orbital_Couplings_with_Energy_Hot.png）
couplings = pd.DataFrame()
for i in range(1, 8):
    for j in range(i + 1, 8):
        couplings[f'f{i}*f{j}'] = data[f'f{i}'] * data[f'f{j}']
couplings['Energy'] = data['Energy']

# 计算耦合特征与能量的相关性
coupling_correlation = couplings.corr()['Energy'][:-1].values

# 初始化7x7矩阵为0
coupling_correlation_matrix = pd.DataFrame(0, index=[f'f{i}' for i in range(1, 8)],
                                           columns=[f'f{j}' for j in range(1, 8)])

index = 0
for i in range(1, 8):
    for j in range(i + 1, 8):
        correlation_value = coupling_correlation[index]
        coupling_correlation_matrix.loc[f'f{i}', f'f{j}'] = correlation_value
        coupling_correlation_matrix.loc[f'f{j}', f'f{i}'] = correlation_value
        index += 1

# 绘制热力图
plt.figure(figsize=(10, 8))
sns.heatmap(coupling_correlation_matrix.astype(float), annot=True, cmap='coolwarm', fmt=".2f", annot_kws={"size": 10})
plt.title('Correlation of f-Orbital Couplings with Energy')
plt.xlabel('f-Orbital')
plt.ylabel('f-Orbital')
plt.show()

# 主成分分析（PCA）
features = ['f1', 'f2', 'f3', 'f4', 'f5', 'f6', 'f7']
x = data[features]
y = data['Energy']

# 标准化
scaler = StandardScaler()
x_scaled = scaler.fit_transform(x)

# PCA降维
pca = PCA(n_components=3)
principal_components = pca.fit_transform(x_scaled)
pca_df = pd.DataFrame(data=principal_components, columns=['PC1', 'PC2', 'PC3'])
pca_df['Energy'] = y.values

# 自定义颜色映射
cmap = mcolors.LinearSegmentedColormap.from_list("custom", ["red", "yellow", "green", "blue"])

# 绘制PCA结果3D散点图
fig = plt.figure(figsize=(16, 10))
ax = fig.add_subplot(111, projection='3d')

# 绘制散点
sc = ax.scatter(pca_df['PC1'], pca_df['PC2'], pca_df['PC3'], c=pca_df['Energy'], cmap=cmap)

# 添加颜色条
cbar = plt.colorbar(sc, ax=ax, pad=0.1)
cbar.set_label('Energy')

# 设置轴标签
ax.set_xlabel('PC1')
ax.set_ylabel('PC2')
ax.set_zlabel('PC3')
ax.set_title('PCA of f-Orbital Occupations (3D)')

plt.show()