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

plt.rcParams['font.family'] = 'SimHei'  # 设置中文字体

# 创建数据集
data = {
    '标号': range(1, 32),
    'diameter': [8.3, 8.6, 8.8, 10.5, 10.7, 10.8, 11.0, 11.0, 11.1, 11.2, 11.3, 11.4, 11.4, 11.7, 
                 12.0, 12.9, 12.9, 13.3, 13.7, 13.8, 14.0, 14.2, 14.5, 16.0, 16.3, 17.3, 17.5, 17.9, 
                 18.0, 18.0, 20.6],
    'height': [70, 65, 63, 72, 81, 83, 66, 75, 80, 75, 79, 76, 76, 69, 75, 74, 85, 86, 71, 64, 78, 
               80, 74, 72, 77, 81, 82, 80, 80, 80, 87],
    'volume': [10.3, 10.3, 10.2, 16.4, 18.8, 19.7, 15.6, 18.2, 22.6, 19.9, 24.2, 21.0, 21.4, 21.3, 
               19.1, 22.2, 33.8, 27.4, 25.7, 24.9, 34.5, 31.7, 36.3, 38.3, 42.6, 55.4, 55.7, 58.3, 
               51.5, 51.0, 77.0]
}

# 先创建 DataFrame
df = pd.DataFrame(data)
df['diameter'] = df['diameter'] ** 2
p = 3
n = 31
e = np.ones((31, 1))

# 提取自变量矩阵 X（添加常数项）
X = np.column_stack((np.ones((31, 1)), df['diameter'], df['height']))

# 提取因变量 Y
Y = df['volume'].values.reshape(-1, 1)

# 最小二乘解 β = (XᵀX)^(-1) XᵀY
Xt = X.T
XtX = Xt @ X
XtX_inv = np.linalg.inv(XtX)

H = X @ XtX_inv @ Xt
I = np.eye(H.shape[0])

beta = XtX_inv @ Xt @ Y
print("回归系数 β：")
print(beta)

# Y预测值
Y_cal = X @ beta

# SSE和sigema计算
epsilon = Y - Y_cal
SSE = epsilon.T @ (I - H) @ epsilon
print("SSE：", SSE)
sigema = SSE / (n-p)
print("误差方差：", sigema)
MSE = SSE / (n - p)
print("MSE:", MSE)

# SST
SST = Y.T @ (I- (e @ e.T) / n) @ Y
print("SST：", SST)


# SSR
SSR = Y.T @ H @ (I- (e @ e.T) / n) @ H @ Y
print("SSR:", SSR)
MSR = SSR / (p - 1)
print("MSR:", MSR)


# R2计算
R2 = 1 - SSE / SST
print("复相关系数的平方：", R2)