未记录事实核验 · 状态：维护中

记录仅说明核验时间与维护状态，不保证适用于所有环境；使用前请留意正文中的版本、适用条件与风险提示。

📈 线性回归（Linear Regression）

🤔 什么是线性回归？

生活例子：

• 根据房屋面积预测房价：面积越大，房价越高
• 根据学习时间预测考试成绩：学习时间越长，成绩越好
• 根据广告投入预测销售额：广告费越多，销售越好

📐 线性回归的数学原理

简单线性回归

模型方程：$y = \beta_0 + \beta_1 x + \varepsilon$

其中：

• $y$：因变量（响应变量）
• $x$：自变量（解释变量）
• $\beta_0$：截距（当x=0时的线性预测值；若E(ε|X)=0，则为对应的条件均值）
• $\beta_1$：斜率（x每增加1个单位，y的变化量）
• $\varepsilon$：误差项

参数估计（最小二乘法）

斜率：$\hat{\beta_1} = \frac{\sum_{i=1}^{n}(x_i - \bar{x})(y_i - \bar{y})}{\sum_{i=1}^{n}(x_i - \bar{x})^2}$

截距：$\hat{\beta_0} = \bar{y} - \hat{\beta_1}\bar{x}$

import numpy as np
import matplotlib.pyplot as plt

def simple_linear_regression(x, y):
    """手动实现简单线性回归"""
    n = len(x)
    x_mean = np.mean(x)
    y_mean = np.mean(y)
    
 # 计算斜率和截距
    numerator = np.sum((x - x_mean) * (y - y_mean))
    denominator = np.sum((x - x_mean) ** 2)
    beta_1 = numerator / denominator
    beta_0 = y_mean - beta_1 * x_mean
    
 # 预测值
    y_pred = beta_0 + beta_1 * x
    
 # 计算R²
    ss_res = np.sum((y - y_pred) ** 2)  # 残差平方和
    ss_tot = np.sum((y - y_mean) ** 2)  # 总平方和
    r_squared = 1 - (ss_res / ss_tot)
    
    return beta_0, beta_1, y_pred, r_squared

# 示例数据：学习时间 vs 考试成绩
np.random.seed(42)
study_hours = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
exam_scores = 50 + 4 * study_hours + np.random.normal(0, 3, 10)

# 拟合模型
intercept, slope, predictions, r2 = simple_linear_regression(study_hours, exam_scores)

print(f"回归方程: 成绩 = {intercept:.2f} + {slope:.2f} × 学习时间")
print(f"R² = {r2:.3f}")
print(f"解释: 在本模拟样本的线性拟合中，学习时间每增加1小时，拟合成绩平均增加{slope:.2f}分；不据此推断因果效果")

# 可视化
plt.figure(figsize=(10, 6))
plt.scatter(study_hours, exam_scores, color='blue', label='实际数据')
plt.plot(study_hours, predictions, color='red', label=f'回归线 (R²={r2:.3f})')
plt.xlabel('学习时间 (小时)')
plt.ylabel('考试成绩 (分)')
plt.title('学习时间 vs 考试成绩')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

💻 Python实现

1. 使用scikit-learn

from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score, mean_squared_error
import numpy as np
import pandas as pd

# 创建数据
np.random.seed(42)
X = np.random.uniform(1, 10, 100).reshape(-1, 1)  # 自变量需要是2D数组
y = 2 * X.flatten() + 3 + np.random.normal(0, 1, 100)  # 真实关系: y = 2x + 3 + 噪声

# 创建和训练模型
model = LinearRegression()
model.fit(X, y)

# 预测
y_pred = model.predict(X)

# 模型评估
r2 = r2_score(y, y_pred)
mse = mean_squared_error(y, y_pred)
rmse = np.sqrt(mse)

print(f"模型参数:")
print(f"截距 (β₀): {model.intercept_:.3f}")
print(f"斜率 (β₁): {model.coef_[0]:.3f}")
print(f"\n模型评估:")
print(f"R²: {r2:.3f}")
print(f"RMSE: {rmse:.3f}")

# 预测新值
new_X = np.array([[5], [7], [9]])
new_pred = model.predict(new_X)
print(f"\n预测结果:")
for i, (x_val, pred) in enumerate(zip(new_X.flatten(), new_pred)):
    print(f"X = {x_val}: 预测y = {pred:.2f}")

2. 使用statsmodels（统计分析）

import statsmodels.api as sm
import pandas as pd

# 准备数据
df = pd.DataFrame({'X': X.flatten(), 'y': y})

# 添加常数项（截距）
X_with_const = sm.add_constant(df['X'])

# 拟合模型
model_stats = sm.OLS(df['y'], X_with_const).fit()

# 详细统计报告
print(model_stats.summary())

# 获取关键统计量
print(f"\n关键统计量:")
print(f"R²: {model_stats.rsquared:.3f}")
print(f"调整R²: {model_stats.rsquared_adj:.3f}")
print(f"F统计量: {model_stats.fvalue:.3f}")
print(f"P值: {model_stats.f_pvalue:.6f}")

# 参数的置信区间
conf_int = model_stats.conf_int()
print(f"\n参数95%置信区间:")
print(conf_int)

📊 模型评估指标

1. 决定系数（R²）

def calculate_r_squared(y_true, y_pred):
    """计算R²的详细步骤"""
    y_mean = np.mean(y_true)
    
 # 总平方和 (Total Sum of Squares)
    ss_tot = np.sum((y_true - y_mean) ** 2)
    
 # 残差平方和 (Residual Sum of Squares)
    ss_res = np.sum((y_true - y_pred) ** 2)
    
 # 回归平方和 (Explained Sum of Squares)
    ss_reg = ss_tot - ss_res
    
 # R²
    r_squared = 1 - (ss_res / ss_tot)
 # 或者 r_squared = ss_reg / ss_tot
    
    print(f"总变异 (SS_tot): {ss_tot:.2f}")
    print(f"残差变异 (SS_res): {ss_res:.2f}")
    print(f"解释变异 (SS_reg): {ss_reg:.2f}")
    print(f"R² = 1 - SS_res/SS_tot = {r_squared:.3f}")
    print(f"解释: 模型解释了{r_squared*100:.1f}%的y变异")
    
    return r_squared

# 使用之前的数据
r2_manual = calculate_r_squared(y, y_pred)

2. 其他评估指标

from sklearn.metrics import mean_absolute_error, mean_squared_error, mean_absolute_percentage_error

def regression_metrics(y_true, y_pred):
    """计算各种回归评估指标"""
    mae = mean_absolute_error(y_true, y_pred)
    mse = mean_squared_error(y_true, y_pred)
    rmse = np.sqrt(mse)
    mape = mean_absolute_percentage_error(y_true, y_pred)
    r2 = r2_score(y_true, y_pred)
    
    print(f"回归评估指标:")
    print(f"MAE (平均绝对误差): {mae:.3f}")
    print(f"MSE (均方误差): {mse:.3f}")
    print(f"RMSE (均方根误差): {rmse:.3f}")
    print(f"MAPE (平均绝对百分比误差): {mape:.2%}")
    print(f"R² (决定系数): {r2:.3f}")
    
    return {'MAE': mae, 'MSE': mse, 'RMSE': rmse, 'MAPE': mape, 'R2': r2}

# 评估模型
metrics = regression_metrics(y, y_pred)

🔍 模型诊断

1. 残差分析

def residual_analysis(y_true, y_pred, X):
    """残差分析"""
    residuals = y_true - y_pred
    
    fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(12, 10))
    
 # 1. 残差 vs 拟合值
    ax1.scatter(y_pred, residuals, alpha=0.7)
    ax1.axhline(y=0, color='red', linestyle='--')
    ax1.set_xlabel('拟合值')
    ax1.set_ylabel('残差')
    ax1.set_title('残差 vs 拟合值')
    
 # 2. Q-Q图（检验正态性）
    from scipy import stats
    stats.probplot(residuals, dist="norm", plot=ax2)
    ax2.set_title('残差Q-Q图')
    
 # 3. 残差直方图
    ax3.hist(residuals, bins=20, density=True, alpha=0.7)
    ax3.set_xlabel('残差')
    ax3.set_ylabel('密度')
    ax3.set_title('残差分布')
    
 # 4. 残差 vs 自变量
    ax4.scatter(X.flatten(), residuals, alpha=0.7)
    ax4.axhline(y=0, color='red', linestyle='--')
    ax4.set_xlabel('X值')
    ax4.set_ylabel('残差')
    ax4.set_title('残差 vs X值')
    
    plt.tight_layout()
    plt.show()
    
 # 统计检验
 # Shapiro-Wilk正态性检验
    shapiro_stat, shapiro_p = stats.shapiro(residuals)
    print(f"Shapiro-Wilk正态性检验:")
    print(f"统计量: {shapiro_stat:.4f}, p值: {shapiro_p:.4f}")
    if shapiro_p > 0.05:
        print("在0.05显著性水平下未拒绝正态性原假设，不能据此证明残差服从正态分布")
    else:
        print("在0.05显著性水平下拒绝正态性原假设，存在偏离正态分布的统计证据")

# 进行残差分析
residual_analysis(y, y_pred, X)

2. 异常值检测

def detect_regression_outliers(X, y, model):
    """检测回归中的异常值"""
    y_pred = model.predict(X)
    residuals = y - y_pred
    
 # 标准化残差
    residual_std = np.sqrt(np.sum(residuals ** 2) / (len(y) - X.shape[1] - 1))
    standardized_residuals = residuals / residual_std
    
 # Cook's Distance (影响点)
    n = len(y)
    p = X.shape[1] + 1  # 参数个数（包括截距）
    
 # 杠杆值 (Leverage)
    X_with_const = np.column_stack([np.ones(len(X)), X])
    H = X_with_const @ np.linalg.inv(X_with_const.T @ X_with_const) @ X_with_const.T
    leverage = np.diag(H)
    
 # Cook's Distance
    cooks_d = (standardized_residuals ** 2) * (leverage / (1 - leverage) ** 2) / p
    
 # 异常值标准
    outlier_threshold = 2  # 标准化残差阈值
    leverage_threshold = 2 * p / n  # 杠杆值阈值
    cooks_threshold = 4 / n  # Cook's Distance阈值
    
    outliers = []
    for i in range(len(y)):
        reasons = []
        if abs(standardized_residuals[i]) > outlier_threshold:
            reasons.append(f"高残差({standardized_residuals[i]:.2f})")
        if leverage[i] > leverage_threshold:
            reasons.append(f"高杠杆({leverage[i]:.3f})")
        if cooks_d[i] > cooks_threshold:
            reasons.append(f"高影响({cooks_d[i]:.3f})")
        
        if reasons:
            outliers.append((i, X[i].flatten()[0], y[i], reasons))
    
    print(f"异常值检测结果:")
    print(f"标准化残差阈值: ±{outlier_threshold}")
    print(f"杠杆值阈值: {leverage_threshold:.3f}")
    print(f"Cook's Distance阈值: {cooks_threshold:.3f}")
    print(f"\n发现{len(outliers)}个潜在异常值:")
    
    for idx, x_val, y_val, reasons in outliers:
        print(f"第{idx}个观测: X={x_val:.2f}, y={y_val:.2f}, 原因: {', '.join(reasons)}")
    
    return outliers, standardized_residuals, leverage, cooks_d

# 检测异常值
outliers, std_res, lev, cooks = detect_regression_outliers(X, y, model)

🏠 实际应用案例

案例1：房价预测

# 生成房价数据
np.random.seed(42)
house_size = np.random.uniform(50, 200, 100)  # 房屋面积 (平方米)
# 价格与面积成正比，加上一些噪声
house_price = 10000 + 300 * house_size + np.random.normal(0, 5000, 100)

# 数据预处理
from sklearn.model_selection import train_test_split

X_house = house_size.reshape(-1, 1)
y_house = house_price

# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(
    X_house, y_house, test_size=0.2, random_state=42
)

# 训练模型
house_model = LinearRegression()
house_model.fit(X_train, y_train)

# 预测
y_train_pred = house_model.predict(X_train)
y_test_pred = house_model.predict(X_test)

# 评估
train_r2 = r2_score(y_train, y_train_pred)
test_r2 = r2_score(y_test, y_test_pred)

print(f"房价预测模型:")
print(f"模型方程: 价格 = {house_model.intercept_:.0f} + {house_model.coef_[0]:.0f} × 面积")
print(f"训练集R²: {train_r2:.3f}")
print(f"测试集R²: {test_r2:.3f}")

# 可视化
plt.figure(figsize=(12, 5))

plt.subplot(1, 2, 1)
plt.scatter(X_train, y_train, alpha=0.6, label='训练数据')
plt.scatter(X_test, y_test, alpha=0.6, label='测试数据', color='orange')
x_range = np.linspace(X_house.min(), X_house.max(), 100).reshape(-1, 1)
y_range = house_model.predict(x_range)
plt.plot(x_range, y_range, 'r-', label='回归线')
plt.xlabel('房屋面积 (㎡)')
plt.ylabel('房价 (元)')
plt.title('房屋面积 vs 房价')
plt.legend()

plt.subplot(1, 2, 2)
plt.scatter(y_test, y_test_pred, alpha=0.6)
plt.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], 'r--')
plt.xlabel('实际房价')
plt.ylabel('预测房价')
plt.title(f'预测效果 (R²={test_r2:.3f})')

plt.tight_layout()
plt.show()

# 实际预测
print(f"\n房价预测示例:")
test_sizes = [80, 120, 150]
for size in test_sizes:
    predicted_price = house_model.predict(np.array([[size]]))[0]
    print(f"{size}㎡房屋预测价格: {predicted_price:,.0f}元")

案例2：销售预测

# 生成销售数据
np.random.seed(42)
ad_spend = np.random.uniform(10, 100, 80)  # 广告投入 (万元)
# 模拟销售额与广告投入相关，正二次项使条件均值的边际效应递增
sales = 50 + 2 * ad_spend + 0.01 * ad_spend**2 + np.random.normal(0, 10, 80)

# 线性回归分析
X_sales = ad_spend.reshape(-1, 1)
y_sales = sales

sales_model = LinearRegression()
sales_model.fit(X_sales, y_sales)
y_sales_pred = sales_model.predict(X_sales)

# 模型评估
sales_r2 = r2_score(y_sales, y_sales_pred)
sales_rmse = np.sqrt(mean_squared_error(y_sales, y_sales_pred))

print(f"销售预测模型:")
print(f"模型方程: 销售额 = {sales_model.intercept_:.1f} + {sales_model.coef_[0]:.2f} × 广告投入")
print(f"R²: {sales_r2:.3f}")
print(f"RMSE: {sales_rmse:.1f}万元")

# ROI分析
marginal_return = sales_model.coef_[0]  # 边际收益
print(f"\n投资回报分析:")
print(f"在本模拟样本的线性拟合中，广告投入每增加1万元，拟合销售额平均增加{marginal_return:.2f}万元；这不是因果增量收益估计")
print(f"仅假设拟合斜率代表因果增量、增量销售额全部转为收益且除广告费外无其他成本时，简化边际ROI为{(marginal_return - 1) * 100:.1f}%；本示例未验证这些条件")

# 可视化
plt.figure(figsize=(10, 6))
plt.scatter(ad_spend, sales, alpha=0.6, label='实际数据')
x_range = np.linspace(ad_spend.min(), ad_spend.max(), 100).reshape(-1, 1)
y_range = sales_model.predict(x_range)
plt.plot(x_range, y_range, 'r-', label=f'线性回归 (R²={sales_r2:.3f})')
plt.xlabel('广告投入 (万元)')
plt.ylabel('销售额 (万元)')
plt.title('广告投入 vs 销售额')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

⚠️ 线性回归的假设和限制

基本假设

1. 线性关系：X和Y之间存在线性关系
2. 独立性：观测值相互独立
3. 同方差性：误差项的方差恒定
4. 误差条件均值为零与正态性边界：解释回归条件均值通常要求E(ε|X)=0；经典有限样本t/F推断还依赖误差正态性等条件，OLS拟合本身不要求误差服从正态分布
5. 无多重共线性：自变量间无强相关

违反假设的处理

# 生成非线性数据
np.random.seed(42)
x_nonlinear = np.linspace(0, 10, 100)
y_nonlinear = 2 * x_nonlinear**2 + np.random.normal(0, 10, 100)

# 普通线性回归
X_nl = x_nonlinear.reshape(-1, 1)
linear_model = LinearRegression()
linear_model.fit(X_nl, y_nonlinear)
y_linear_pred = linear_model.predict(X_nl)

# 多项式回归
from sklearn.preprocessing import PolynomialFeatures
poly_features = PolynomialFeatures(degree=2)
X_poly = poly_features.fit_transform(X_nl)
poly_model = LinearRegression()
poly_model.fit(X_poly, y_nonlinear)
y_poly_pred = poly_model.predict(X_poly)

# 比较效果
linear_r2 = r2_score(y_nonlinear, y_linear_pred)
poly_r2 = r2_score(y_nonlinear, y_poly_pred)

plt.figure(figsize=(12, 5))

plt.subplot(1, 2, 1)
plt.scatter(x_nonlinear, y_nonlinear, alpha=0.6, label='实际数据')
plt.plot(x_nonlinear, y_linear_pred, 'r-', label=f'线性回归 (R²={linear_r2:.3f})')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('线性回归处理非线性数据')
plt.legend()

plt.subplot(1, 2, 2)
plt.scatter(x_nonlinear, y_nonlinear, alpha=0.6, label='实际数据')
plt.plot(x_nonlinear, y_poly_pred, 'g-', label=f'多项式回归 (R²={poly_r2:.3f})')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('多项式回归处理非线性数据')
plt.legend()

plt.tight_layout()
plt.show()

print(f"线性回归R²: {linear_r2:.3f}")
print(f"多项式回归R²: {poly_r2:.3f}")
print(f"改进: {(poly_r2 - linear_r2):.3f}")

🛠️ 实用技巧

1. 特征工程

# 创建多特征数据
np.random.seed(42)
n_samples = 200

# 原始特征
feature1 = np.random.uniform(1, 10, n_samples)
feature2 = np.random.uniform(20, 50, n_samples)
feature3 = np.random.uniform(100, 500, n_samples)

# 目标变量
y_multi = (2 * feature1 + 
           0.5 * feature2 + 
           0.01 * feature3 + 
           0.1 * feature1 * feature2 +  # 交互项
           np.random.normal(0, 5, n_samples))

# 原始特征矩阵
X_original = np.column_stack([feature1, feature2, feature3])

# 特征工程
from sklearn.preprocessing import StandardScaler, PolynomialFeatures

# 1. 标准化
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X_original)

# 2. 添加交互项和多项式特征
poly = PolynomialFeatures(degree=2, interaction_only=True, include_bias=False)
X_poly = poly.fit_transform(X_original)

# 3. 比较不同特征组合的效果
models = {
    '原始特征': LinearRegression().fit(X_original, y_multi),
    '标准化特征': LinearRegression().fit(X_scaled, y_multi),
    '多项式特征': LinearRegression().fit(X_poly, y_multi)
}

print("不同特征工程方法的效果:")
for name, model in models.items():
    if name == '原始特征':
        y_pred = model.predict(X_original)
    elif name == '标准化特征':
        y_pred = model.predict(X_scaled)
    else:
        y_pred = model.predict(X_poly)
    
    r2 = r2_score(y_multi, y_pred)
    print(f"{name}: R² = {r2:.3f}")

2. 交叉验证

from sklearn.model_selection import cross_val_score, KFold

# 使用交叉验证评估模型稳定性
cv_scores = cross_val_score(LinearRegression(), X_original, y_multi, 
                           cv=KFold(n_splits=5, shuffle=True, random_state=42),
                           scoring='r2')

print(f"5折交叉验证结果:")
print(f"R²分数: {cv_scores}")
print(f"平均R²: {cv_scores.mean():.3f} ± {cv_scores.std():.3f}")

# 学习曲线
from sklearn.model_selection import learning_curve

train_sizes, train_scores, val_scores = learning_curve(
    LinearRegression(), X_original, y_multi, 
    train_sizes=np.linspace(0.1, 1.0, 10),
    cv=5, scoring='r2'
)

plt.figure(figsize=(10, 6))
plt.plot(train_sizes, train_scores.mean(axis=1), 'o-', label='训练分数')
plt.plot(train_sizes, val_scores.mean(axis=1), 'o-', label='验证分数')
plt.fill_between(train_sizes, train_scores.mean(axis=1) - train_scores.std(axis=1),
                 train_scores.mean(axis=1) + train_scores.std(axis=1), alpha=0.2)
plt.fill_between(train_sizes, val_scores.mean(axis=1) - val_scores.std(axis=1),
                 val_scores.mean(axis=1) + val_scores.std(axis=1), alpha=0.2)
plt.xlabel('训练样本数')
plt.ylabel('R²分数')
plt.title('学习曲线')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

🔗 相关知识

• 相关
• 皮尔逊相关系数
• Presto/连接器
• 多元线性回归
• 逻辑回归
• 假设检验

#统计学 #线性回归 #回归分析 #预测 #机器学习 #基础