import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)
# 4 种典型关系
x = np.linspace(0, 10, 50)
fig, axes = plt.subplots(1, 4, figsize=(16, 4))
for ax, slope, noise, title in zip(
    axes, [2, -1.5, 0, 0],
    [1, 1, 0.5, 3],
    ["正相关", "负相关", "无线性相关", "曲线关系"]):
    y = slope * x + np.random.normal(0, noise, len(x))
    ax.scatter(x, y, alpha=0.6)
    ax.set_title(title)
    ax.grid(True, alpha=0.3)
plt.tight_layout(); plt.show()

from scipy.stats import pearsonr
np.random.seed(42)
x = np.random.normal(0, 1, 100)
y = 0.7 * x + np.random.normal(0, 0.5, 100)
r, p = pearsonr(x, y)
print(f"r = {r:.3f}, p = {p:.4f}")
# r ≈ 0.81 (强), p ≈ 0 (极显著)

# scipy 的 pearsonr 已经给出 p 值
r, p = pearsonr(x, y)  # p 就是 ρ=0 的检验 p 值

import numpy as np
from sklearn.linear_model import LinearRegression

# 数据
np.random.seed(42)
x = np.random.uniform(0, 10, 50)
y = 3 + 2 * x + np.random.normal(0, 2, 50)

# sklearn
model = LinearRegression()
model.fit(x.reshape(-1, 1), y)
print(f"b0 = {model.intercept_:.3f}, b1 = {model.coef_[0]:.3f}")
# b0 ≈ 3, b1 ≈ 2 (跟真实 3 + 2x 吻合!)

import matplotlib.pyplot as plt

y_pred = model.predict(x.reshape(-1, 1))
residuals = y - y_pred

fig, axes = plt.subplots(1, 3, figsize=(15, 4))
# 1. 残差 vs 拟合值 (应随机散布)
axes[0].scatter(y_pred, residuals, alpha=0.6)
axes[0].axhline(0, color="red", linestyle="--")
axes[0].set_xlabel("拟合值"); axes[0].set_ylabel("残差")
axes[0].set_title("残差 vs 拟合值")

# 2. 残差 Q-Q 图 (应近似直线)
from scipy.stats import probplot
probplot(residuals, dist="norm", plot=axes[1])
axes[1].set_title("Q-Q 图 (正态性)")

# 3. 残差直方图
axes[2].hist(residuals, bins=15, edgecolor="black")
axes[2].set_title("残差分布")
plt.tight_layout(); plt.show()

# 数值检验
from scipy.stats import shapiro, bartlett
print(f"Shapiro 残差正态性: p = {shapiro(residuals).pvalue:.3f}")
# p > 0.05 → 残差正态, 假设成立

print(f"R² = {model.score(x.reshape(-1, 1), y):.3f}")
# 0.85 左右, 拟合好

from scipy.stats import t

# 计算预测区间
n = len(x)
x_mean = x.mean()
sse = np.sum(residuals**2)
s_e = np.sqrt(sse / (n - 2))  # 残差标准误
t_crit = t.ppf(0.975, df=n-2)

x_new = np.array([[5.0]])
y_pred = model.predict(x_new)[0]
# 预测区间
pi_se = s_e * np.sqrt(1 + 1/n + (5 - x_mean)**2 / np.sum((x - x_mean)**2))
pi = (y_pred - t_crit * pi_se, y_pred + t_crit * pi_se)
print(f"x=5 时预测值 = {y_pred:.2f}, 95% PI = [{pi[0]:.2f}, {pi[1]:.2f}]")

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt

# 模拟房价数据
np.random.seed(42)
n = 200
size = np.random.uniform(50, 200, n)  # 面积 m²
rooms = (size / 30 + np.random.normal(0, 0.5, n)).astype(int)  # 房间数
age = np.random.uniform(0, 30, n)       # 房龄
price = (size * 0.5 + rooms * 5 - age * 0.3 +
         np.random.normal(0, 10, n))    # 万元

df = pd.DataFrame({"size": size, "rooms": rooms, "age": age, "price": price})

# 1. 探索性分析
print("相关系数:")
print(df.corr()["price"].round(3))

# 2. 拟合简单回归
X = df[["size"]]
y = df["price"]
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

model = LinearRegression().fit(X_train, y_train)
y_pred = model.predict(X_test)
print(f"\n单变量: R² = {r2_score(y_test, y_pred):.3f}, RMSE = {np.sqrt(mean_squared_error(y_test, y_pred)):.2f}")

# 3. 多元回归
model2 = LinearRegression().fit(df[["size", "rooms", "age"]], y)
y_pred2 = model2.predict(df[["size", "rooms", "age"]])
print(f"多变量: R² = {r2_score(y, y_pred2):.3f}")

# 4. 系数解释
print("\n回归系数:")
for name, coef in zip(["size", "rooms", "age"], model2.coef_):
    print(f"  {name:6s}: {coef:+.3f}")
# size: +0.50 (每 +1 m², 房价 +0.5 万)
# rooms: +5.0 (每多 1 房间, +5 万)
# age: -0.3 (每老 1 年, -0.3 万)

import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score
from scipy import stats

def full_regression_diagnosis(X, y, feature_names=None):
    """完整回归诊断: 拟合 + 残差 + 异常值 + 显著检验"""
    if feature_names is None:
        feature_names = [f"X{i}" for i in range(X.shape[1])]
    
    model = LinearRegression().fit(X, y)
    y_pred = model.predict(X)
    residuals = y - y_pred
    
    print("="*60)
    print(f"R² = {r2_score(y, y_pred):.4f}, 调整 R² = {1 - (1-r2_score(y, y_pred))*(len(y)-1)/(len(y)-X.shape[1]-1):.4f}")
    print(f"系数: {dict(zip(feature_names, model.coef_.round(4)))}")
    print(f"截距: {model.intercept_:.4f}")
    print("="*60)
    
    # 残差诊断
    print(f"\n残差诊断:")
    print(f"  均值 = {residuals.mean():.4f} (应≈0)")
    print(f"  标准差 = {residuals.std(ddof=1):.4f}")
    print(f"  Shapiro 正态性 p = {stats.shapiro(residuals).pvalue:.4f}")
    
    # 异常值检测
    infl = model.predict(X)  # 占位
    student_resid = residuals / residuals.std(ddof=1)
    outliers = np.abs(student_resid) > 3
    print(f"  异常值 (|z| > 3): {outliers.sum()} 个")
    
    # 系数显著检验 (近似 t 检验)
    n, k = X.shape
    mse = np.sum(residuals**2) / (n - k - 1)
    X_with_intercept = np.column_stack([np.ones(n), X])
    var_coef = mse * np.linalg.inv(X_with_intercept.T @ X_with_intercept).diagonal()
    se = np.sqrt(var_coef)
    t_stats = np.concatenate([[model.intercept_], model.coef_]) / se
    p_vals = 2 * (1 - stats.t.cdf(np.abs(t_stats), df=n-k-1))
    
    print(f"\n系数显著检验:")
    print(f"  {'变量':15s} {'系数':>10s} {'SE':>10s} {'t':>10s} {'p':>10s}")
    for name, b, s, t_v, p_v in zip(
        ["截距"] + list(feature_names),
        np.concatenate([[model.intercept_], model.coef_]),
        se, t_stats, p_vals):
        sig = "***" if p_v < 0.001 else ("**" if p_v < 0.01 else ("*" if p_v < 0.05 else ""))
        print(f"  {name:15s} {b:10.4f} {s:10.4f} {t_v:10.3f} {p_v:10.4f} {sig}")
    
    return model

# 测试
np.random.seed(42)
X = np.column_stack([
    np.random.normal(0, 1, 100),
    np.random.normal(0, 1, 100),
    np.random.normal(0, 1, 100),
])
y = 1 + 0.5 * X[:, 0] + 0.3 * X[:, 1] + np.random.normal(0, 0.5, 100)

full_regression_diagnosis(X, y, ["X1", "X2", "X3"])