第24章 线性回归

import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler, PolynomialFeatures
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.metrics import r2_score, mean_squared_error

# 生成模拟房价数据
np.random.seed(42)
n = 200

# 原始特征
area = np.random.uniform(50, 200, n)        # 面积
age = np.random.uniform(0, 50, n)           # 房龄
rooms = np.random.randint(1, 6, n)          # 房间数
distance = np.random.uniform(1, 30, n)      # 距市中心距离

# 真实关系包含：
# 1. 面积的平方效应（递减）
# 2. 面积×房间数交互
# 3. 距离的对数效应
price = (100 + 
         3*area - 0.005*area**2 +        # 面积非线性
         10*rooms +                       # 房间数
         0.05*area*rooms -                # 交互项
         5*np.log(distance+1) -           # 距离对数
         0.5*age +                        # 房龄
         np.random.randn(n)*10)           # 噪声

# 创建DataFrame
df = pd.DataFrame({
    '面积': area,
    '房龄': age,
    '房间数': rooms,
    '距离': distance,
    '价格': price
})

print("数据前5行:")
print(df.head())
print(f"\n数据形状: {df.shape}")

# 特征工程
def feature_engineering(df, fit=True, scaler=None, poly=None):
    """完整的特征工程流程"""
    X = df[['面积', '房龄', '房间数', '距离']].copy()
    
    # 1. 对数变换（距离）
    X['距离_log'] = np.log(X['距离'] + 1)
    
    # 2. 多项式特征（面积二次项）
    X['面积平方'] = X['面积'] ** 2
    
    # 3. 交互项（面积×房间数）
    X['面积×房间数'] = X['面积'] * X['房间数']
    
    # 移除原始距离（用对数版本替代）
    X = X.drop('距离', axis=1)
    
    # 4. 标准化
    if fit:
        scaler = StandardScaler()
        X_scaled = scaler.fit_transform(X)
    else:
        X_scaled = scaler.transform(X)
    
    return X_scaled, scaler

# 划分数据
train_df, test_df = train_test_split(df, test_size=0.2, random_state=42)
y_train = train_df['价格'].values
y_test = test_df['价格'].values

# 应用特征工程
X_train, scaler = feature_engineering(train_df, fit=True)
X_test, _ = feature_engineering(test_df, fit=False, scaler=scaler)

# 对比模型
print("\n" + "="*70)
print("模型对比")
print("="*70)

# 基线：只用原始特征
X_train_base = train_df[['面积', '房龄', '房间数', '距离']].values
X_test_base = test_df[['面积', '房龄', '房间数', '距离']].values
scaler_base = StandardScaler()
X_train_base = scaler_base.fit_transform(X_train_base)
X_test_base = scaler_base.transform(X_test_base)

model_base = LinearRegression()
model_base.fit(X_train_base, y_train)
y_pred_base = model_base.predict(X_test_base)
r2_base = r2_score(y_test, y_pred_base)
rmse_base = np.sqrt(mean_squared_error(y_test, y_pred_base))

print(f"基线模型（原始特征）:")
print(f"  测试集 R² = {r2_base:.4f}")
print(f"  测试集 RMSE = {rmse_base:.2f}")

# 特征工程模型
model_fe = LinearRegression()
model_fe.fit(X_train, y_train)
y_pred_fe = model_fe.predict(X_test)
r2_fe = r2_score(y_test, y_pred_fe)
rmse_fe = np.sqrt(mean_squared_error(y_test, y_pred_fe))

print(f"\n特征工程模型:")
print(f"  测试集 R² = {r2_fe:.4f}")
print(f"  测试集 RMSE = {rmse_fe:.2f}")
print(f"\n改进:")
print(f"  R²提升: {(r2_fe - r2_base):.4f}")
print(f"  RMSE降低: {(rmse_base - rmse_fe):.2f}")

# 特征重要性
feature_names = ['面积', '房龄', '房间数', '距离_log', '面积平方', '面积×房间数']
importances = np.abs(model_fe.coef_)
indices = np.argsort(importances)[::-1]

print(f"\n特征重要性（按系数绝对值）:")
for i in indices:
    print(f"  {feature_names[i]:<15}: {importances[i]:.3f}")

print("="*70)

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import statsmodels.api as sm
from sklearn.linear_model import LinearRegression

# 生成数据
np.random.seed(42)
n = 100
X = np.random.uniform(0, 10, n)
true_beta0, true_beta1 = 5, 2
y = true_beta0 + true_beta1 * X + np.random.randn(n) * 3

# 方法1：使用statsmodels（推荐，提供完整统计信息）
X_with_const = sm.add_constant(X)  # 添加截距列
model_sm = sm.OLS(y, X_with_const)
results = model_sm.fit()

print("="*70)
print("statsmodels完整输出")
print("="*70)
print(results.summary())
print("\n" + "="*70)

# 提取关键统计量
print("\n关键统计量：")
print("-"*70)
print(f"系数估计:")
print(f"  β₀ (截距): {results.params[0]:.4f}")
print(f"  β₁ (斜率): {results.params[1]:.4f}")
print(f"\n标准误:")
print(f"  SE(β₀): {results.bse[0]:.4f}")
print(f"  SE(β₁): {results.bse[1]:.4f}")
print(f"\nt统计量和p值:")
print(f"  β₀: t={results.tvalues[0]:.4f}, p={results.pvalues[0]:.4e}")
print(f"  β₁: t={results.tvalues[1]:.4f}, p={results.pvalues[1]:.4e}")
print(f"\n95%置信区间:")
conf_int = results.conf_int(alpha=0.05)
print(f"  β₀: [{conf_int[0][0]:.4f}, {conf_int[0][1]:.4f}]")
print(f"  β₁: [{conf_int[1][0]:.4f}, {conf_int[1][1]:.4f}]")
print(f"\nR²和F统计量:")
print(f"  R² = {results.rsquared:.4f}")
print(f"  调整R² = {results.rsquared_adj:.4f}")
print(f"  F统计量 = {results.fvalue:.4f}, p值 = {results.f_pvalue:.4e}")
print(f"\n残差标准差:")
print(f"  σ̂ = {np.sqrt(results.mse_resid):.4f}")
print("="*70)

# 预测与预测区间
print("\n" + "="*70)
print("预测与预测区间")
print("="*70)

# 新样本点
X_new = np.array([2, 5, 8])
X_new_const = sm.add_constant(X_new)

# 点预测
y_pred = results.predict(X_new_const)

# 置信区间（均值的区间）
pred_ci = results.get_prediction(X_new_const).conf_int(alpha=0.05)

# 预测区间（个体的区间）
pred_pi = results.get_prediction(X_new_const).conf_int(obs=True, alpha=0.05)

print(f"\n{'X':<10} {'预测值':<12} {'95% CI':<30} {'95% PI':<30}")
print("-"*70)
for i, x in enumerate(X_new):
    print(f"{x:<10.1f} {y_pred[i]:<12.2f} "
          f"[{pred_ci[i,0]:.2f}, {pred_ci[i,1]:.2f}]{'':>12} "
          f"[{pred_pi[i,0]:.2f}, {pred_pi[i,1]:.2f}]")
print("="*70)

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 1. 拟合曲线与置信带
ax1 = axes[0, 0]
X_plot = np.linspace(X.min(), X.max(), 200)
X_plot_const = sm.add_constant(X_plot)
y_plot = results.predict(X_plot_const)
pred_plot = results.get_prediction(X_plot_const)
ci_plot = pred_plot.conf_int(alpha=0.05)
pi_plot = pred_plot.conf_int(obs=True, alpha=0.05)

ax1.scatter(X, y, alpha=0.5, s=30, label='数据')
ax1.plot(X_plot, y_plot, 'r-', linewidth=2, label='拟合线')
ax1.fill_between(X_plot, ci_plot[:, 0], ci_plot[:, 1], 
                  alpha=0.3, color='blue', label='95% 置信带（均值）')
ax1.fill_between(X_plot, pi_plot[:, 0], pi_plot[:, 1], 
                  alpha=0.2, color='green', label='95% 预测带（个体）')
ax1.set_xlabel('X')
ax1.set_ylabel('y')
ax1.set_title('拟合曲线、置信带与预测带')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 2. 系数的置信区间
ax2 = axes[0, 1]
params = results.params
conf_int = results.conf_int(alpha=0.05)
param_names = ['β₀\n(截距)', 'β₁\n(斜率)']
x_pos = np.arange(len(params))

ax2.errorbar(x_pos, params, 
             yerr=[params - conf_int[:, 0], conf_int[:, 1] - params],
             fmt='o', markersize=10, capsize=10, capthick=2, linewidth=2)
ax2.axhline(y=0, color='red', linestyle='--', linewidth=2, label='H₀: β=0')
ax2.set_xticks(x_pos)
ax2.set_xticklabels(param_names)
ax2.set_ylabel('系数值')
ax2.set_title('系数的95%置信区间')
ax2.legend()
ax2.grid(True, alpha=0.3, axis='y')

# 3. t统计量
ax3 = axes[1, 0]
t_values = results.tvalues
colors = ['green' if abs(t) > 2 else 'red' for t in t_values]
ax3.bar(param_names, t_values, color=colors, alpha=0.7)
ax3.axhline(y=2, color='black', linestyle='--', linewidth=1, label='临界值≈2')
ax3.axhline(y=-2, color='black', linestyle='--', linewidth=1)
ax3.axhline(y=0, color='gray', linestyle='-', linewidth=1)
ax3.set_ylabel('t统计量')
ax3.set_title('t统计量（|t|>2→显著）')
ax3.legend()
ax3.grid(True, alpha=0.3, axis='y')

# 4. 残差图
ax4 = axes[1, 1]
residuals = results.resid
fitted = results.fittedvalues
ax4.scatter(fitted, residuals, alpha=0.5, s=30)
ax4.axhline(y=0, color='red', linestyle='--', linewidth=2)
ax4.set_xlabel('拟合值')
ax4.set_ylabel('残差')
ax4.set_title('残差图')
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# 方法2：手动计算（教学用）
print("\n" + "="*70)
print("手动计算（验证）")
print("="*70)

# 使用sklearn拟合
from sklearn.linear_model import LinearRegression
model_sk = LinearRegression()
model_sk.fit(X.reshape(-1, 1), y)

# 预测和残差
y_pred_train = model_sk.predict(X.reshape(-1, 1))
residuals = y - y_pred_train
RSS = np.sum(residuals**2)
TSS = np.sum((y - y.mean())**2)

# 自由度
n = len(y)
p = 1  # 1个自变量
df_resid = n - p - 1

# 残差标准差
sigma_hat = np.sqrt(RSS / df_resid)

# 标准误（简化计算，假设X已中心化）
X_centered = X - X.mean()
var_x = np.sum(X_centered**2)
SE_beta1 = sigma_hat / np.sqrt(var_x)
SE_beta0 = sigma_hat * np.sqrt(1/n + X.mean()**2/var_x)

# t统计量
t_beta0 = model_sk.intercept_ / SE_beta0
t_beta1 = model_sk.coef_[0] / SE_beta1

print(f"\n手动计算结果:")
print(f"  σ̂ = {sigma_hat:.4f}")
print(f"  SE(β₀) = {SE_beta0:.4f}")
print(f"  SE(β₁) = {SE_beta1:.4f}")
print(f"  t(β₀) = {t_beta0:.4f}")
print(f"  t(β₁) = {t_beta1:.4f}")
print(f"\n与statsmodels对比：基本一致✓")
print("="*70)