import numpy as np
import pandas as pd
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA

# 从图像中读取数据
# 创建一个包含49位女性血糖数据的DataFrame
# X1, X2, X3 表示空腹情况下三个不同时刻的血糖含量
# Y1, Y2, Y3 表示摄入等量食糖一小时后三个时刻的血糖含量

# 第一部分数据（1-16号）
data1 = [
    [60, 69, 62, 97, 68, 98],
    [56, 53, 84, 103, 78, 107],
    [80, 66, 76, 66, 99, 130],
    [55, 58, 70, 90, 85, 114],
    [62, 75, 68, 116, 130, 91],
    [74, 84, 70, 96, 101, 103],
    [64, 71, 66, 77, 102, 130],
    [73, 70, 73, 82, 110, 109],
    [68, 67, 72, 76, 85, 115],
    [69, 82, 61, 96, 133, 127],
    [60, 67, 61, 150, 134, 121],
    [70, 74, 76, 150, 158, 100],
    [66, 74, 78, 150, 131, 142],
    [83, 76, 74, 162, 98, 105],
    [68, 66, 90, 119, 85, 109],
    [78, 63, 75, 123, 98, 138]
]

# 第二部分数据（17-32号）
data2 = [
    [103, 77, 77, 160, 117, 121],
    [77, 68, 74, 144, 71, 153],
    [66, 77, 68, 77, 82, 91],
    [70, 70, 72, 114, 93, 122],
    [72, 85, 71, 77, 87, 109],
    [91, 74, 93, 118, 115, 130],
    [89, 75, 73, 170, 147, 121],
    [75, 82, 76, 153, 132, 115],
    [84, 71, 66, 143, 107, 100],
    [76, 70, 64, 114, 113, 129],
    [75, 90, 86, 75, 94, 116],
    [74, 77, 80, 116, 81, 97],
    [89, 71, 69, 63, 87, 70],
    [78, 75, 80, 105, 132, 80],
    [74, 66, 71, 85, 113, 113],
    [71, 80, 76, 81, 87, 86]
]

# 第三部分数据（33-48号）
data3 = [
    [69, 75, 73, 120, 89, 95],
    [90, 103, 74, 107, 109, 101],
    [90, 76, 81, 99, 111, 98],
    [48, 77, 75, 113, 124, 97],
    [66, 93, 97, 136, 142, 122],
    [74, 70, 78, 109, 88, 105],
    [70, 74, 71, 72, 90, 75],
    [65, 75, 66, 130, 101, 90],
    [66, 80, 86, 130, 132, 144],
    [77, 67, 74, 83, 92, 107],
    [70, 67, 100, 150, 142, 146],
    [73, 76, 81, 119, 120, 119],
    [78, 90, 77, 122, 145, 149],
    [73, 68, 80, 102, 99, 122],
    [72, 83, 77, 104, 69, 96]
]

# 第四部分数据（49号）
data4 = [
    [52, 70, 76, 92, 94, 100]
]

# 合并所有数据
data = data1 + data2 + data3 + data4

# 创建DataFrame
columns = ['X1', 'X2', 'X3', 'Y1', 'Y2', 'Y3']
df = pd.DataFrame(data, columns=columns)

print("原始数据:")
print(df.head())
print(f"数据形状: {df.shape}")

# 题目要求:
# 1. 从样本协方差矩阵S出发作主成分分析
print("\n1. 基于样本协方差矩阵S的主成分分析:")

# 标准化数据
scaler = StandardScaler()
scaled_data = scaler.fit_transform(df)

# 计算协方差矩阵S (默认使用协方差矩阵)
pca_S = PCA(n_components=6)
pca_S.fit(scaled_data)

# 特征值
eigenvalues_S = pca_S.explained_variance_
print("\n特征值 (基于S):")
for i, val in enumerate(eigenvalues_S):
    print(f"λ{i+1} = {val:.4f}")

# 主成分贡献率
contribution_ratios_S = pca_S.explained_variance_ratio_
cumulative_ratios_S = np.cumsum(contribution_ratios_S)

print("\n各主成分贡献率 (基于S):")
for i, ratio in enumerate(contribution_ratios_S):
    print(f"主成分 {i+1}: {ratio:.4f} ({ratio*100:.2f}%)")

print("\n累计贡献率 (基于S):")
for i, ratio in enumerate(cumulative_ratios_S):
    print(f"前 {i+1} 个主成分: {ratio:.4f} ({ratio*100:.2f}%)")

# 主成分系数 (特征向量)
components_S = pca_S.components_
print("\n主成分系数 (基于S):")
component_df_S = pd.DataFrame(components_S, 
                            index=[f'PC{i+1}' for i in range(len(components_S))],
                            columns=columns)
print(component_df_S)

# 计算主成分得分
pc_scores_S = pca_S.transform(scaled_data)
print("\n前三个主成分得分 (部分示例, 基于S):")
pc_scores_df_S = pd.DataFrame(pc_scores_S[:, :3], 
                            columns=['PC1', 'PC2', 'PC3'],
                            index=range(1, len(data)+1))
print(pc_scores_df_S.head())

# 2. 从样本相关系数矩阵R出发作主成分分析
print("\n2. 基于样本相关系数矩阵R的主成分分析:")

# 计算相关系数矩阵R
correlation_matrix = df.corr()
print("\n相关系数矩阵R:")
print(correlation_matrix)

# 使用相关系数矩阵进行PCA (设置svd_solver='full'确保使用相关系数矩阵)
pca_R = PCA(n_components=6, svd_solver='full')
# 此时scaled_data已经是标准化后的数据，协方差矩阵等同于相关系数矩阵
pca_R.fit(scaled_data)

# 特征值
eigenvalues_R = pca_R.explained_variance_
print("\n特征值 (基于R):")
for i, val in enumerate(eigenvalues_R):
    print(f"λ{i+1} = {val:.4f}")

# 主成分贡献率
contribution_ratios_R = pca_R.explained_variance_ratio_
cumulative_ratios_R = np.cumsum(contribution_ratios_R)

print("\n各主成分贡献率 (基于R):")
for i, ratio in enumerate(contribution_ratios_R):
    print(f"主成分 {i+1}: {ratio:.4f} ({ratio*100:.2f}%)")

print("\n累计贡献率 (基于R):")
for i, ratio in enumerate(cumulative_ratios_R):
    print(f"前 {i+1} 个主成分: {ratio:.4f} ({ratio*100:.2f}%)")

# 主成分系数 (特征向量)
components_R = pca_R.components_
print("\n主成分系数 (基于R):")
component_df_R = pd.DataFrame(components_R, 
                            index=[f'PC{i+1}' for i in range(len(components_R))],
                            columns=columns)
print(component_df_R)

# 计算主成分得分
pc_scores_R = pca_R.transform(scaled_data)
print("\n前三个主成分得分 (部分示例, 基于R):")
pc_scores_df_R = pd.DataFrame(pc_scores_R[:, :3], 
                            columns=['PC1', 'PC2', 'PC3'],
                            index=range(1, len(data)+1))
print(pc_scores_df_R.head())

# 结果分析和解释
print("\n主成分保留分析:")

# 分析S矩阵的结果
print("\n基于S矩阵:")
# 计算保留至少85%信息需要的主成分数
for i, ratio in enumerate(cumulative_ratios_S):
    if ratio >= 0.85:
        print(f"保留前 {i+1} 个主成分可解释 {ratio*100:.2f}% 的总方差")
        break

# 分析R矩阵的结果
print("\n基于R矩阵:")
# 计算保留至少85%信息需要的主成分数
for i, ratio in enumerate(cumulative_ratios_R):
    if ratio >= 0.85:
        print(f"保留前 {i+1} 个主成分可解释 {ratio*100:.2f}% 的总方差")
        break

# 分析前两个主成分的含义
print("\n主成分解释:")
print("基于R矩阵的前两个主成分解释:")
for i, comp in enumerate([components_R[0], components_R[1]]):
    print(f"主成分 {i+1} 的系数:")
    for j, col in enumerate(columns):
        print(f"  {col}: {comp[j]:.4f}")
    
    # 简单解释
    if i == 0:
        print("  解释: 第一主成分可能代表总体血糖水平")
    else:
        print("  解释: 第二主成分可能表示空腹与进食后血糖变化的对比")

# 比较S和R的结果，判断哪个更合理
print("\n分析比较:")
print("基于协方差矩阵S和相关系数矩阵R的主成分分析结果的比较:")
print("1. 当原始变量的量纲/方差差异较大时，基于R的分析通常更合理，因为它消除了量纲影响。")
print("2. 对于本例中的血糖数据，各指标单位相同(mg/100ml)，但空腹和进食后的血糖含量水平和波动范围可能有较大差异。")
print("3. 如果各变量的方差差异较大，则基于S的分析可能会偏向方差大的变量，而基于R的分析则给予各变量同等重要性。")
print("4. 通常在多数应用中，基于R的主成分分析更为常用，因为它避免了因变量量纲不同或方差差异过大导致的偏倚。")
print("5. 对于此题，由于血糖数据的取值范围较广，空腹和餐后血糖的方差可能存在差异，因此基于R的主成分分析可能更合理。")