import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.formula.api import ols
import matplotlib.pyplot as plt



# 1. 设计参数
num_factors = 3
num_levels = 3

# 2. 手动定义 L9(3^4) 正交表
# L9 表是固定的，表示 9 次试验，4 列（因素），每列 3 个水平。 可以做3因素3水平设计，1列是空白列，也可以做4因素3水平。
# 我们只需要其中的 3 列（L, M, N）
# 水平值使用 1, 2, 3
orthogonal_design_l9 = np.array([
    [1, 1, 1], # 试验 1: L1 M1 N1
    [1, 2, 2], # 试验 2: L1 M2 N2
    [1, 3, 3], # 试验 3: L1 M3 N3
    [2, 1, 2], # 试验 4: L2 M1 N2
    [2, 2, 3], # 试验 5: L2 M2 N3
    [2, 3, 1], # 试验 6: L2 M3 N1
    [3, 1, 3], # 试验 7: L3 M1 N3
    [3, 2, 1], # 试验 8: L3 M2 N1
    [3, 3, 2]  # 试验 9: L3 M3 N2
])

# 3. 可视化设计
print("--- L9 正交试验设计 (3因素 3水平) ---")
print("因素: L, M, N (水平值: 1, 2, 3)")
print(orthogonal_design_l9)
print(f"\n试验次数: {orthogonal_design_l9.shape[0]}")




# 4. 正交试验分析（数据收集和分析框架）
# 假设的实验结果 (响应值 Y)
Y = np.array([9.5, 12.1, 9.8, 15.0, 16.5, 12.9, 10.0, 13.2, 10.8])

# 将设计矩阵和结果合并
df = pd.DataFrame(orthogonal_design_l9, columns=['L', 'M', 'N'])
df['Y'] = Y

print("\n\n\n--- 原始数据表 ---")
print(df)



# --- 极差分析 (Range Analysis) 框架 ---

def range_analysis(df, factors):
    """进行简单的极差分析"""
    results = {}
    
    for factor in factors:
        # 按因素和水平分组，计算响应值 Y 的平均值 (K 值)
        K_values = df.groupby(factor)['Y'].mean()
        
        # 计算极差 R = Max(K) - Min(K)
        R = K_values.max() - K_values.min()
        
        results[factor] = {
            'K_values': K_values.to_dict(),
            'Range (R)': R
        }
        
    # 按极差 R 排序，R 越大，影响越大
    ranked_factors = sorted(results.keys(), key=lambda k: results[k]['Range (R)'], reverse=True)
    
    print("\n\n\n--- 极差分析结果 ---")
    for factor in ranked_factors:
        print(f"因素 {factor}:")
        # 格式化输出 K 值，保留两位小数
        k_str = {k: f"{v:.2f}" for k, v in results[factor]['K_values'].items()}
        print(f"  各水平平均值 (K): {k_str}")
        print(f"  极差 (R): {results[factor]['Range (R)']:.2f}")

    print(f"\n影响显著性排序 (R值): {' > '.join(ranked_factors)}")
    
    # 确定最佳水平组合 (最大化 Y)
    best_combination = {}
    for factor in factors:
        K_values = results[factor]['K_values']
        # 找到 K 值最大的水平
        best_level = max(K_values, key=K_values.get)
        best_combination[factor] = int(best_level)
        
    print(f"\n最佳水平组合 (最大化 Y): {best_combination}")
    
# 执行极差分析
range_analysis(df, factors=['L', 'M', 'N'])




# --- 方差分析 (ANOVA) 框架 ---


print("\n\n\n--- 方差分析 (ANOVA) 框架 ---")

# 将因素 A, B, C 转换为 categorical 类型以便进行方差分析
df_anova = df.copy()
df_anova['L'] = df_anova['L'].astype('category')
df_anova['M'] = df_anova['M'].astype('category')
df_anova['N'] = df_anova['N'].astype('category')

# 构建模型公式: Y ~ C(L) + C(M) + C(N) 表示 Y 是 L, M, N 的函数
model = ols('Y ~ C(L) + C(M) + C(N)', data=df_anova).fit()

# 输出 ANOVA 表
# 注意：L9 表的自由度限制（总共 8 个自由度）使得无法包含所有交互项。
# 此处模型将所有未被 L, M, N 主效应解释的变异都归入 Residual（残差/误差）。
anova_table = sm.stats.anova_lm(model, typ=1)

print(anova_table)




# --- 绘制主效应图 ---

mean_effects = {}
for factor in ['L', 'M', 'N']:
    mean_effects[factor] = df.groupby(factor)['Y'].mean()

print("\n\n\n--- 主效应（平均响应）---")
for factor, means in mean_effects.items():
    print(f"{factor}: {means}")


fig, axes = plt.subplots(1, 3, figsize=(12, 4))
for i, factor in enumerate(['L', 'M', 'N']):
    mean_effects[factor].plot(ax=axes[i], marker='o')
    axes[i].set_title(f'Main Effect of {factor}')
    axes[i].set_xlabel('Level')
    axes[i].set_ylabel('Mean Response')
plt.tight_layout()
plt.show()