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

# 创建数据
workers = ['W1', 'W2', 'W3', 'W4']
machines = ['A1', 'A2', 'A3']

# 产量数据
production = np.array([
    [50, 47, 47, 53],  # A1
    [63, 54, 57, 58],  # A2
    [52, 42, 41, 48]   # A3
])

# 创建DataFrame以便于分析
data = []
for i, machine in enumerate(machines):
    for j, worker in enumerate(workers):
        data.append({
            'Machine': machine,
            'Worker': worker,
            'Production': production[i, j]
        })

df = pd.DataFrame(data)

# 1. 建立方差分析表
print("两因素方差分析：")
model = ols('Production ~ C(Machine) + C(Worker)', data=df).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
print(anova_table)

# 判断显著性
alpha = 0.05
print(f"\n在显著性水平 α = {alpha} 下的检验结果：")
print(f"机床因素的p值: {anova_table.loc['C(Machine)', 'PR(>F)']:.4f}")
if anova_table.loc['C(Machine)', 'PR(>F)'] < alpha:
    print("结论：不同机床的日产量有显著差异")
else:
    print("结论：不同机床的日产量无显著差异")

print(f"工人因素的p值: {anova_table.loc['C(Worker)', 'PR(>F)']:.4f}")
if anova_table.loc['C(Worker)', 'PR(>F)'] < alpha:
    print("结论：不同工人的日产量有显著差异")
else:
    print("结论：不同工人的日产量无显著差异")

# 2. 计算各因素水平的平均值
worker_means = df.groupby('Worker')['Production'].mean()
machine_means = df.groupby('Machine')['Production'].mean()

print("\n各工人的平均日产量：")
print(worker_means)
print("\n各机床的平均日产量：")
print(machine_means)

# Bonferroni同时置信区间计算函数
def bonferroni_intervals(group_means, mse, n_per_group, alpha=0.05, num_comparisons=None):
    """
    计算Bonferroni同时置信区间
    
    参数:
    group_means: 各组均值
    mse: 均方误差
    n_per_group: 每组的样本量
    alpha: 显著性水平
    num_comparisons: 总比较次数，如果为None则自动计算
    
    返回:
    包含各组间比较结果的列表
    """
    groups = group_means.index
    k = len(groups)
    
    if num_comparisons is None:
        num_comparisons = k * (k - 1) // 2
    
    # 残差自由度 = 总样本数 - 组数 (对于单因素设计)
    # 对于本题的两因素设计，我们从ANOVA表中获取残差自由度
    dof = anova_table.loc['Residual', 'df']
    
    # Bonferroni调整的t临界值
    t_critical = stats.t.ppf(1 - alpha/(2*num_comparisons), dof)
    
    # 组间差异的标准误
    se = np.sqrt(2 * mse / n_per_group)
    
    intervals = []
    for i in range(k):
        for j in range(i+1, k):
            diff = group_means.iloc[i] - group_means.iloc[j]
            margin = t_critical * se
            lower = diff - margin
            upper = diff + margin
            
            # 检查是否显著
            significant = (lower > 0) or (upper < 0)
            
            intervals.append({
                'group1': groups[i],
                'group2': groups[j],
                'diff': diff,
                'lower': lower,
                'upper': upper,
                'significant': significant
            })
    
    return intervals

# 计算总比较次数
total_comparisons = (len(workers) * (len(workers) - 1) // 2) + (len(machines) * (len(machines) - 1) // 2)

# 计算MSE（均方误差）
mse = anova_table.loc['Residual', 'sum_sq'] / anova_table.loc['Residual', 'df']

# 题目要求的是置信度不小于90%的Bonferroni同时置信区间，对应alpha=0.10
alpha_confidence = 0.10

# 工人Bonferroni同时置信区间 (每个工人对应3个机床数据)
worker_intervals = bonferroni_intervals(worker_means, mse, n_per_group=3, 
                                         alpha=alpha_confidence, num_comparisons=total_comparisons)

print("\n工人平均日产量的Bonferroni 90%同时置信区间：")
for interval in worker_intervals:
    print(f"{interval['group1']} vs {interval['group2']}: 差值 = {interval['diff']:.2f}, "
          f"置信区间 = [{interval['lower']:.2f}, {interval['upper']:.2f}]")
    print(f"  结论: {'有显著差异' if interval['significant'] else '无显著差异'}")

# 机床Bonferroni同时置信区间 (每台机床对应4个工人数据)
machine_intervals = bonferroni_intervals(machine_means, mse, n_per_group=4, 
                                          alpha=alpha_confidence, num_comparisons=total_comparisons)

print("\n机床平均日产量的Bonferroni 90%同时置信区间：")
for interval in machine_intervals:
    print(f"{interval['group1']} vs {interval['group2']}: 差值 = {interval['diff']:.2f}, "
          f"置信区间 = [{interval['lower']:.2f}, {interval['upper']:.2f}]")
    print(f"  结论: {'有显著差异' if interval['significant'] else '无显著差异'}")

# 可视化结果
plt.figure(figsize=(12, 10))

# 工人日产量箱线图
plt.subplot(2, 2, 1)
df.boxplot(column='Production', by='Worker')
plt.title('各工人的日产量分布')
plt.ylabel('日产量(件)')
plt.xlabel('工人')
plt.suptitle('')

# 机床日产量箱线图
plt.subplot(2, 2, 2)
df.boxplot(column='Production', by='Machine')
plt.title('各机床的日产量分布')
plt.ylabel('日产量(件)')
plt.xlabel('机床')
plt.suptitle('')

# 工人均值条形图
plt.subplot(2, 2, 3)
x = np.arange(len(workers))
plt.bar(x, worker_means.values, width=0.6, color='skyblue')
# 计算每个工人的置信区间宽度
worker_ci_width = np.zeros(len(workers))
for interval in worker_intervals:
    if interval['group1'] == 'W1' and interval['group2'] == 'W2':
        # 使用第一个区间的宽度作为估计
        margin = (interval['upper'] - interval['lower']) / 2
        worker_ci_width = np.ones(len(workers)) * margin

plt.errorbar(x, worker_means.values, yerr=worker_ci_width, fmt='none', color='black', capsize=5)
plt.xticks(x, workers)
plt.title('工人平均日产量及90%置信区间')
plt.ylabel('平均日产量(件)')
plt.xlabel('工人')

# 机床均值条形图
plt.subplot(2, 2, 4)
x = np.arange(len(machines))
plt.bar(x, machine_means.values, width=0.6, color='lightgreen')
# 计算每个机床的置信区间宽度
machine_ci_width = np.zeros(len(machines))
for interval in machine_intervals:
    if interval['group1'] == 'A1' and interval['group2'] == 'A2':
        # 使用第一个区间的宽度作为估计
        margin = (interval['upper'] - interval['lower']) / 2
        machine_ci_width = np.ones(len(machines)) * margin

plt.errorbar(x, machine_means.values, yerr=machine_ci_width, fmt='none', color='black', capsize=5)
plt.xticks(x, machines)
plt.title('机床平均日产量及90%置信区间')
plt.ylabel('平均日产量(件)')
plt.xlabel('机床')

plt.tight_layout()
plt.show()

# 分析结果总结
print("\n综合分析结果：")
print("1. 在α=0.05的显著性水平下，不同机床的日产量存在显著差异(p值=0.0008)")
print("2. 在α=0.05的显著性水平下，不同工人的日产量存在显著差异(p值=0.0220)")
print("3. Bonferroni 90%同时置信区间分析表明：")
print("   - 工人方面：W1与W2的日产量存在显著差异，其他工人之间无显著差异")
print("   - 机床方面：A2与A1、A2与A3的日产量存在显著差异，A1与A3之间无显著差异")
print("4. 机床A2的平均日产量最高(58.00件)，机床A3的平均日产量最低(45.75件)")
print("5. 工人W1的平均日产量最高(55.00件)，工人W2的平均日产量最低(47.67件)")