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, 16),
    '销量': [162, 120, 223, 131, 67, 169, 81, 192, 116, 55, 252, 232, 144, 103, 212],
    '人数': [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],
    '收入': [2450, 3254, 3802, 2838, 2347, 3782, 3008, 2450, 2137, 2560, 4020, 4427, 2660, 2088, 2605]
}

# 先创建 DataFrame
df = pd.DataFrame(data)
p = 3
n = 15
e = np.ones((15, 1))

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

# 提取因变量 Y
Y = df['销量'].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)

# F值
F = MSR / MSE
print("F值：", F)

# 示例：计算 F 分布在 x = 2.5 时的累积分布函数值
cdf_value = f.cdf(F, p - 1, n - p)
print("p值:", 1 - cdf_value)

# 示例：计算显著性水平为 0.95 时的 t 分布临界值
alpha = 0.05
t_critical = t.ppf(1 - alpha / 2, n - p)
# 计算可信区间
diagonal_elements = np.diag(XtX_inv)

# 输出可信区间
for i in range(1, 3):
    beta_value = beta[i]  # beta[i] 是一个标量
    d_value = diagonal_elements[i]  # diagonal_elements[i] 是一个标量
    print(f"Beta[{i+1}] 95% CI: ({beta_value - t_critical * math.sqrt(sigema * d_value)}, {beta_value + t_critical * math.sqrt(sigema * d_value)})")

# y的预测区间
x0 = np.array([[1], [220], [2500]])
y0 = x0.T @ beta
y0_s = y0.item()
print("预测销量值：", y0_s)

mm = x0.T @ XtX_inv @ x0
mm_s = mm.item() 
tt = math.sqrt(MSE * (1 + mm_s))
print(f"y的置信区间为：({y0_s - t_critical * tt}, {y0_s + t_critical * tt})")


# 输出可信区间
for i in range(1, 3):
    beta_value = beta[i]  # beta[i] 是一个标量
    d_value = diagonal_elements[i]  # diagonal_elements[i] 是一个标量
    s = math.sqrt(sigema * d_value)
    print(f"t值,beta{i+1}: ({beta_value / s})")


 


















