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

# 设置中文显示
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False    # 用来正常显示负号

# 数据输入
low = [7.6, 8.2, 6.8, 5.8, 6.9, 6.6, 6.3, 7.7, 6.0]
medium = [6.7, 8.1, 9.4, 8.6, 7.8, 7.7, 8.9, 7.9, 8.3, 8.7, 7.1, 8.4]
high = [8.5, 9.7, 10.1, 7.8, 9.6, 9.5]

# 创建数据框
data = {'投入水平': ['低'] * len(low) + ['中'] * len(medium) + ['高'] * len(high),
        '生产能力提高量': low + medium + high}
df = pd.DataFrame(data)

# 描述性统计
print("描述性统计：")
desc_stats = df.groupby('投入水平')['生产能力提高量'].describe()
print(desc_stats)

# (1) 建立方差分析表
print("\n(1) 方差分析表:")
# 使用statsmodels进行方差分析
formula = '生产能力提高量 ~ C(投入水平)'
model = ols(formula, data=df).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
print(anova_table)

# F检验结果分析
alpha = 0.05
p_value = anova_table['PR(>F)'][0]
print(f"\n在显著性水平α={alpha}下:")
if p_value < alpha:
    print(f"P值({p_value:.4f}) < {alpha}，拒绝原假设，不同科研经费投入对当年生产力的提高有显著影响。")
else:
    print(f"P值({p_value:.4f}) >= {alpha}，接受原假设，不同科研经费投入对当年生产力的提高无显著影响。")

# (2) 计算均值和95%置信区间
print("\n(2) 各组均值及95%置信区间:")
group_means = df.groupby('投入水平')['生产能力提高量'].mean()
mu_L = group_means['低']
mu_M = group_means['中']
mu_H = group_means['高']

# 计算各组的95%置信区间
n_L = len(low)
n_M = len(medium)
n_H = len(high)

# 计算标准差
std_L = np.std(low, ddof=1)
std_M = np.std(medium, ddof=1)
std_H = np.std(high, ddof=1)

# 计算标准误
se_L = std_L / np.sqrt(n_L)
se_M = std_M / np.sqrt(n_M)
se_H = std_H / np.sqrt(n_H)

# t临界值
t_L = stats.t.ppf(1 - alpha/2, n_L - 1)
t_M = stats.t.ppf(1 - alpha/2, n_M - 1)
t_H = stats.t.ppf(1 - alpha/2, n_H - 1)

# 计算置信区间
CI_L = (mu_L - t_L * se_L, mu_L + t_L * se_L)
CI_M = (mu_M - t_M * se_M, mu_M + t_M * se_M)
CI_H = (mu_H - t_H * se_H, mu_H + t_H * se_H)

print(f"低投入组均值μ_L = {mu_L:.4f}, 95%置信区间: [{CI_L[0]:.4f}, {CI_L[1]:.4f}]")
print(f"中投入组均值μ_M = {mu_M:.4f}, 95%置信区间: [{CI_M[0]:.4f}, {CI_M[1]:.4f}]")
print(f"高投入组均值μ_H = {mu_H:.4f}, 95%置信区间: [{CI_H[0]:.4f}, {CI_H[1]:.4f}]")

# 计算Bonferroni同时置信区间
# 对于三个均值，有3个对比：μ_L-μ_M, μ_L-μ_H, μ_M-μ_H
# Bonferroni调整后的α为α/3
alpha_adj = alpha / 3

# MSE (均方误差)从方差分析表获取
MSE = anova_table['sum_sq'][1] / anova_table['df'][1]

# 计算差值及其置信区间
diff_LM = mu_L - mu_M
diff_LH = mu_L - mu_H
diff_MH = mu_M - mu_H

# 计算Bonferroni调整后的t临界值
df_error = anova_table['df'][1]  # 误差自由度
t_bonf = stats.t.ppf(1-alpha_adj/2, df_error)

# 计算差值的标准误
se_diff_LM = np.sqrt(MSE * (1/n_L + 1/n_M))
se_diff_LH = np.sqrt(MSE * (1/n_L + 1/n_H))
se_diff_MH = np.sqrt(MSE * (1/n_M + 1/n_H))

# 计算差值的置信区间
CI_diff_LM = (diff_LM - t_bonf * se_diff_LM, diff_LM + t_bonf * se_diff_LM)
CI_diff_LH = (diff_LH - t_bonf * se_diff_LH, diff_LH + t_bonf * se_diff_LH)
CI_diff_MH = (diff_MH - t_bonf * se_diff_MH, diff_MH + t_bonf * se_diff_MH)

print("\nBonferroni同时置信区间 (置信度不小于95%):")
print(f"μ_L - μ_M = {diff_LM:.4f}, 置信区间: [{CI_diff_LM[0]:.4f}, {CI_diff_LM[1]:.4f}]")
print(f"μ_L - μ_H = {diff_LH:.4f}, 置信区间: [{CI_diff_LH[0]:.4f}, {CI_diff_LH[1]:.4f}]")
print(f"μ_M - μ_H = {diff_MH:.4f}, 置信区间: [{CI_diff_MH[0]:.4f}, {CI_diff_MH[1]:.4f}]")

# 判断科研投入对生产能力的影响
print("\n结论分析:")
if p_value < alpha:
    print("根据方差分析结果，不同科研经费投入对生产能力提高量有显著影响。")
    
    # 检查Bonferroni置信区间是否包含0
    significant_pairs = []
    
    if CI_diff_LM[0] > 0 or CI_diff_LM[1] < 0:
        if diff_LM < 0:
            significant_pairs.append("中投入显著高于低投入")
        else:
            significant_pairs.append("低投入显著高于中投入")
    
    if CI_diff_LH[0] > 0 or CI_diff_LH[1] < 0:
        if diff_LH < 0:
            significant_pairs.append("高投入显著高于低投入")
        else:
            significant_pairs.append("低投入显著高于高投入")
    
    if CI_diff_MH[0] > 0 or CI_diff_MH[1] < 0:
        if diff_MH < 0:
            significant_pairs.append("高投入显著高于中投入")
        else:
            significant_pairs.append("中投入显著高于高投入")
    
    if significant_pairs:
        print("具体差异：")
        for pair in significant_pairs:
            print(f"- {pair}")
        
        # 判断是否科研经费投入越高，当年生产能力的改善越显著
        if "高投入显著高于中投入" in significant_pairs and "中投入显著高于低投入" in significant_pairs:
            print("科研经费投入越高，当年生产能力的改善越显著。")
        elif "高投入显著高于低投入" in significant_pairs:
            print("高投入对生产能力提高有显著影响。")
    else:
        print("虽然整体上有显著差异，但在成对比较中未发现显著差异，可能是由于Bonferroni校正过于保守。")
else:
    print("不同科研经费投入对当年生产能力的提高无显著影响。")

# 可视化
plt.figure(figsize=(10, 6))

# 箱线图
plt.subplot(1, 2, 1)
sns.boxplot(x='投入水平', y='生产能力提高量', data=df, order=['低', '中', '高'])
plt.title('不同科研经费投入下生产能力提高量的箱线图')
plt.grid(True, linestyle='--', alpha=0.7)

# 均值及置信区间条形图
plt.subplot(1, 2, 2)
means = [mu_L, mu_M, mu_H]
errors = [(CI_L[1]-CI_L[0])/2, (CI_M[1]-CI_M[0])/2, (CI_H[1]-CI_H[0])/2]
plt.bar(['低', '中', '高'], means, yerr=errors, capsize=10, alpha=0.7)
plt.title('均值及95%置信区间')
plt.ylabel('生产能力提高量')
plt.grid(True, linestyle='--', alpha=0.7)

plt.tight_layout()
plt.show()