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

plt.rcParams['font.family'] = 'Microsoft YaHei'  # 微软雅黑支持负号



# 输入数据
data = {
    "Y": [162, 120, 223, 131, 67, 169, 81, 192, 116, 55, 252, 232, 144, 103, 212],
    "X1": [27.4, 18.0, 37.5, 20.5, 8.6, 26.5, 9.8, 33.0, 19.5, 5.3, 43.0, 37.2, 23.6, 15.7, 37.0],
    "X2": [2450, 3254, 3802, 2838, 2347, 3782, 3008, 2450, 2137, 2560, 4020, 4427, 2660, 2088, 2605]
}

df = pd.DataFrame(data)

# 添加截距项
X = sm.add_constant(df[["X1", "X2"]])
y = df["Y"]

# 拟合模型
model = sm.OLS(y, X).fit()

# 拟合值
fitted_values = model.fittedvalues

# 残差
residuals = model.resid

# 学生化残差
influence = model.get_influence()
studentized_residuals = influence.resid_studentized_internal

# 输出
print("拟合值：\n", fitted_values)
print("残差：\n", residuals)
print("学生化残差：\n", studentized_residuals)

# 绘图设置
plt.figure(figsize=(14, 8))

# 频率直方图
plt.subplot(2, 2, 1)
sns.histplot(studentized_residuals, kde=True)
plt.title("学生化残差的频率分布图")

# QQ图
plt.subplot(2, 2, 2)
sm.qqplot(studentized_residuals, line="45", fit=True, ax=plt.gca())
plt.title("学生化残差的正态QQ图")

# 残差 vs 拟合值图
plt.subplot(2, 2, 3)
plt.scatter(fitted_values, residuals)
plt.axhline(0, color='red', linestyle='--')
plt.xlabel("拟合值")
plt.ylabel("残差")
plt.title("残差 vs 拟合值")

# 学生化残差 vs X1 或 X2（任选一个）
plt.subplot(2, 2, 4)
plt.scatter(df["X1"], studentized_residuals)
plt.axhline(0, color='red', linestyle='--')
plt.xlabel("X1")
plt.ylabel("学生化残差")
plt.title("学生化残差 vs X1")

plt.tight_layout()
plt.show()

# 正态性检验（Shapiro-Wilk）
shapiro_test = stats.shapiro(studentized_residuals)
print("Shapiro-Wilk 检验结果: ", shapiro_test)

# QQ图分位数相关系数
theoretical_quantiles = stats.norm.ppf((np.arange(1, len(studentized_residuals)+1) - 0.5) / len(studentized_residuals))
sorted_residuals = np.sort(studentized_residuals)
correlation = np.corrcoef(theoretical_quantiles, sorted_residuals)[0, 1]
print("分位数相关系数: ", correlation)
