zehua

菜单

Administrator
发布于 2026-09-09 / 0 阅读
0
0

统计学习实验一：线性回归导论

#实验报告#M2 ATSI

实验一：线性回归导论 - 完整详细学习笔记

实验目标与总体设计

这个实验通过三个核心目标系统地介绍了线性回归的理论基础和实践应用。首先是手动实现普通最小二乘法（OLS），深入理解其数学原理和计算过程。其次是通过蒙特卡罗模拟验证Student's t定理，这是线性回归统计推断的理论基础。最后将理论应用到实际数据分析中，展示线性回归在解决真实问题时的考虑因素。

练习1：手动实现OLS

线性回归模型的数据生成机制

线性回归模型假设目标变量 Y 与解释变量 x^{(1)}, ..., x^{(p)} 之间存在线性关系：

Y = \beta_0 + \beta_1 x^{(1)} + \ldots + \beta_p x^{(p)} + \varepsilon

这里 \beta_0 是截距项，代表当所有解释变量为0时的期望响应值；\beta_1, ..., \beta_p 是回归系数，表示各解释变量对响应变量的边际影响；\varepsilon 是随机误差项，满足 \varepsilon \sim \mathcal{N}(0, \sigma^2)，代表模型无法解释的随机变异。

高斯噪声假设是经典线性回归理论的核心。它不仅使得最小二乘估计具有良好的统计性质（BLUE - Best Linear Unbiased Estimator），还为后续的假设检验提供了理论基础。正态分布的选择基于中心极限定理：即使个体误差不服从正态分布，大量独立误差的和也会趋向正态分布。

首先定义模型的均值函数：

def linearModel(X, intercept, coef):
    """
    X: matrix of size p*n WITHOUT THE CONSTANT
    intercept: beta_0
    coef: beta_1...beta_p
    """
    return X @ coef + intercept

这个函数实现了线性模型的核心计算。使用矩阵运算 X @ coef 计算所有样本的线性组合，然后加上截距。这种向量化的实现比循环计算效率高得多，体现了NumPy在科学计算中的优势。

实验参数设置与数据生成

实验采用了最简单的一元线性回归设置：

p = 1           # 一个解释变量
intercept = 1   # β₀ = 1
coef = np.array([2])  # β₁ = 2
sigma = 0.4     # 噪声标准差

选择 p=1 使得模型可以在二维平面上可视化，便于直观理解。真实参数 \beta_0=1, \beta_1=2 意味着真实关系是 Y = 1 + 2X + \varepsilon。噪声标准差 \sigma=0.4 的选择保证了信号与噪声的合理比例，既能看到清晰的线性趋势，又能观察到随机变异的影响。

数据生成过程：

np.random.seed(42)  # 固定随机种子确保可重现性

n = 40
X = np.linspace(-1, 1, n).reshape(-1, 1)  # 在[-1,1]上的规则网格
y = linearModel(X, intercept, coef) + np.random.normal(0, sigma, n)

使用规则网格 np.linspace(-1, 1, n) 而不是随机采样有其特殊考虑。规则网格保证了设计矩阵的良好条件数，避免了数值计算的不稳定性。同时，这种固定设计（fixed design）设置简化了理论分析，是统计学习中的常见假设。reshape(-1, 1) 将一维数组转换为列向量，确保与scikit-learn的接口兼容。

OLS估计的理论基础

普通最小二乘法通过最小化残差平方和来估计参数：

\hat\beta = \arg\min_{\beta} |\underline{y} - \underline{X}\beta|^2

这个优化问题的几何意义是找到使预测值 \underline{X}\beta 与观测值 \underline{y} 之间欧氏距离最小的参数。从投影的角度看，OLS找到了 \underline{y} 在列空间 \text{Col}(\underline{X}) 上的正交投影。

通过对目标函数求导并令其为零，可以得到正规方程（Normal Equation）：

\underline{X}^T\underline{X}\beta = \underline{X}^T\underline{y}

当 \underline{X}^T\underline{X} 可逆时，解析解为：

\hat\beta = (\underline{X}^T\underline{X})^{-1}\underline{X}^T\underline{y}

矩阵 \underline{X}^T\underline{X} 的可逆性要求设计矩阵 \underline{X} 列满秩，即各列线性无关。这在实践中意味着样本量 n 必须大于参数个数 p+1，且不存在完全共线的变量。

OLS的手动实现

实现OLS估计器的完整代码展示了算法的核心步骤：

class OLS:
    def __init__(self, add_intercept=True):
        """
        OLS is an object with the following attributes
        add_intercept: whether to add intercept in the regression
        """
        self.add_intercept = add_intercept  # 是否添加截距
        self.coef = None        # 初始化系数为None（未拟合状态）
        self.intercept = 0      # 初始化截距为0
    
    def fit(self, X, y):
        """
        拟合线性回归模型以求解 β̂ = (X'X)⁻¹X'y
        """
        # 步骤1：构建设计矩阵
        if self.add_intercept:
            # 如果需要截距，在X前面添加一列1
            # np.ones(X.shape[0])创建n个1的向量
            # np.column_stack水平拼接，得到[1, X]的设计矩阵
            X_design = np.column_stack((np.ones(X.shape[0]), X))
        else:
            # 不需要截距，直接使用原始X
            X_design = X
        
        # 步骤2：计算正规方程 β̂ = (X'X)⁻¹X'y
        XtX = X_design.T @ X_design         # 计算 X'X（p×p矩阵）
        Xty = X_design.T @ y                # 计算 X'y（p×1向量）
        beta_hat = np.linalg.inv(XtX) @ Xty # 计算 (X'X)⁻¹X'y
        
        # 步骤3：分离截距和系数
        if self.add_intercept:
            # beta_hat[0]是截距，beta_hat[1:]是其他系数
            self.intercept = beta_hat[0]
            self.coef = beta_hat[1:]
        else:
            # 没有截距，所有都是系数
            self.coef = beta_hat
    
    def predict(self, X):
        """
        对新样本进行预测
        """
        # 检查模型是否已经拟合
        if self.coef is None:
            raise ValueError("Model is not fitted yet.")
        
        # 计算预测值：y = β₀ + Xβ
        y_pred = X @ self.coef
        
        # 如果有截距，加上截距项
        if self.add_intercept:
            y_pred += self.intercept
        
        return y_pred

这个实现有几个关键的设计考虑。首先，add_intercept 参数提供了灵活性，允许用户选择是否包含截距项。在某些特殊情况下（如数据已中心化），可能不需要截距。其次，使用 np.column_stack 构建增广设计矩阵是标准做法，它将常数列添加到原始特征矩阵的第一列。这样做的好处是将截距项统一到参数向量中，简化了矩阵运算。

计算 (X'X)⁻¹ 使用了 np.linalg.inv，这在教学中清晰地展示了公式。但在实际应用中，直接求逆可能存在数值稳定性问题，特别是当 X'X 接近奇异时。更稳健的方法是使用 np.linalg.solve 或QR分解。

模型拟合后，通过比较估计值与真实值验证实现的正确性：

OLSSolver = OLS()
OLSSolver.fit(X, y)
print("OLS estimates:", OLSSolver.intercept, OLSSolver.coef[0])
print("True:", intercept, coef[0])

同时与scikit-learn的实现进行对比，验证结果的一致性：

from sklearn.linear_model import LinearRegression
LRModel = LinearRegression()
LRModel.fit(X, y)
print("OLS with scikit-learn:", LRModel.intercept_, LRModel.coef_[0])

可视化拟合结果

通过绘图直观展示拟合效果，这是评估模型的重要步骤：

fig, ax = plt.subplots(1, 1, figsize=(6, 5))
ax.scatter(X, y, s=8, label='Data')
ax.plot(xplot, yplot, linestyle='dashed', color='red', label='Mean function')
ax.plot(xplot, OLSSolver.predict(xplot), color='black', label='Fitted model')
ax.legend()

使用1000个密集点 xplot = np.linspace(-1, 1, 1000).reshape(-1, 1) 绘制平滑的拟合线。图中展示了三个要素：散点表示观测数据，红色虚线是真实的均值函数，黑色实线是拟合的模型。理想情况下，黑线应该接近红线，偏差反映了估计误差和有限样本的影响。

练习2：Student's t定理的数值验证

Student's t定理的理论内容

Student's t定理是线性回归统计推断的基石，它描述了OLS估计量的抽样分布：

  1. 参数估计的分布：\hat\beta \sim \mathcal{N}\left(\beta, \sigma^2(\underline{X}^T\underline{X})^{-1}\right)

    这表明 \hat\beta 是无偏的（期望值等于真实值），其协方差矩阵为 \sigma^2(\underline{X}^T\underline{X})^{-1}。矩阵 (\underline{X}^T\underline{X})^{-1} 反映了设计矩阵的信息量：\underline{X}^T\underline{X} 越"大"（正定意义下），参数估计的方差越小，估计越精确。

  2. 方差估计的分布：\hat{\sigma}^2 = \frac{|\underline{y}-\underline{X}\hat\beta|^2}{n} \sim \frac{\sigma^2}{n}\chi^2(n-p-1)

    残差平方和服从自由度为 n-p-1 的卡方分布（经过缩放）。自由度 n-p-1 反映了估计 p+1 个参数后剩余的独立信息量。这里使用的是有偏估计（除以 n），无偏估计应除以 n-p-1。

  3. 独立性：\hat\beta 与 \hat{\sigma} 相互独立

    这个性质使得可以构造t统计量进行假设检验。独立性源于正态分布的特殊性质：样本均值与样本方差相互独立。

蒙特卡罗模拟验证

通过2000次重复实验验证理论分布：

# Initialize simulation parameters
n_sim = 2000                    # 模拟次数
beta1 = np.zeros(n_sim)         # 存储每次模拟的β₁估计值
sigma2 = np.zeros(n_sim)        # 存储每次模拟的σ²估计值
OLSSim = OLS()                  # 创建OLS对象用于拟合

# 设置随机种子以确保结果可重现
np.random.seed(42)

# Monte Carlo模拟循环
for i in range(n_sim):
    # 步骤1：生成新的响应变量y
    # linearModel(X, intercept, coef) 生成真实的均值
    # np.random.normal(0, sigma, n) 生成n个独立的噪声项
    y_sim = linearModel(X, intercept, coef) + np.random.normal(0, sigma, n)
    
    # 步骤2：使用OLS拟合模型
    OLSSim.fit(X, y_sim)
    
    # 步骤3：保存β₁的估计值
    # OLSSim.coef[0]是β₁的估计（因为p=1，只有一个系数）
    beta1[i] = OLSSim.coef[0]
    
    # 步骤4：计算并保存σ²的估计值
    # 首先计算拟合值
    y_pred = OLSSim.predict(X)
    
    # 计算残差：实际值 - 预测值
    residuals = y_sim - y_pred
    
    # 计算σ²的估计：残差平方和除以n
    # 注意：这里用的是有偏估计（除以n而不是n-p-1）
    sigma2[i] = np.sum(residuals**2) / n

这个模拟设计体现了固定设计（fixed design）的思想：X 保持不变，只有 y 在每次实验中重新生成。这简化了理论分析，因为 (\underline{X}^T\underline{X})^{-1} 保持不变。每次循环代表一个独立的"实验"，模拟了在相同实验条件下重复收集数据的过程。

理论分布的计算

为了与经验分布比较，需要计算理论分布的参数：

df = n - p - 1  # 自由度
XwithIntercept = np.hstack((np.ones((n, 1)), X))  # 增广设计矩阵
v_beta = np.linalg.inv(XwithIntercept.T @ XwithIntercept)  # (X'X)⁻¹
var_beta1 = sigma**2 * v_beta[1, 1]  # β₁的理论方差

v_beta[1, 1] 提取了对应于 \beta_1 的方差（协方差矩阵的第二个对角元素）。第一个元素 v_beta[0, 0] 对应截距的方差。这个计算展示了设计矩阵如何影响估计精度：如果 X 的变异越大，(\underline{X}^T\underline{X}) 越大，估计方差越小。

分布验证的可视化

通过直方图与理论密度曲线的对比验证Student's t定理：

# β₁的分布
betaplot = coef[0] + 4*np.sqrt(var_beta1) * np.linspace(-1, 1, 1000)
axes[0].hist(beta1, density=True, alpha=0.6, label='Empirical distribution')
axes[0].plot(betaplot, stats.norm.pdf(betaplot, loc=coef, scale=np.sqrt(var_beta1)), 
             color='red', label='Reference')

# σ²的分布
sigplot = np.linspace(0, 1, 1000)
axes[1].hist(sigma2, density=True, alpha=0.6, label='Empirical distribution')
axes[1].plot(sigplot, stats.chi2.pdf(sigplot, df=df, scale=sigma**2/n), 
             color='red', label='Reference')

绘图范围的选择很重要：对于 \hat\beta_1，使用 \pm 4 个标准差覆盖了99.99%的概率质量；对于 \hat\sigma^2，由于卡方分布是非负的，从0开始绘制。density=True 将直方图归一化为概率密度，使其可与理论密度函数比较。

枢轴t统计量

t统计量的构造消除了未知参数 \sigma 的影响，这是进行假设检验的关键：

t_stat = (beta1 - coef[0]) / np.sqrt(n * sigma2 * v_beta[1, 1] / (n - p - 1))

这个统计量的分子是估计误差 \hat\beta_1 - \beta_1，分母是估计的标准误。关键的修正因子 \frac{n}{n-p-1} 将有偏的方差估计转换为无偏估计。在正态假设下，这个统计量精确服从自由度为 n-p-1 的t分布。

t分布的形状介于正态分布和柯西分布之间，具有比正态分布更厚的尾部，反映了用估计方差替代真实方差带来的额外不确定性。当自由度增大时，t分布趋向标准正态分布，这体现了大样本理论。

统计推断的实践工具

实践中通常使用statsmodels包进行统计推断：

import statsmodels.api as sm
XwithIntercept = sm.add_constant(X)
OLSModel = sm.OLS(y, XwithIntercept).fit()
OLSModel.summary()

sm.add_constant(X) 自动在设计矩阵前添加常数列。与scikit-learn不同，statsmodels需要显式添加截距列，这提供了更大的灵活性但需要更多注意。OLSModel.summary() 输出包含了完整的统计信息：参数估计、标准误、t统计量、p值、置信区间、R²、F统计量等。

学生化残差的理论推导

Task 4要求推导学生化残差的分布。残差向量可以表示为：

\hat{\underline{\varepsilon}} = \underline{y} - \underline{X}\hat\beta = (\mathrm{I}_n - P)\underline{\varepsilon}

其中投影矩阵 P = \underline{X}(\underline{X}^T\underline{X})^{-1}\underline{X}^T 将向量投影到 \underline{X} 的列空间。矩阵 \mathrm{I}_n - P 是残差制造矩阵（residual maker matrix），它将向量投影到 \underline{X} 列空间的正交补空间。

投影矩阵的关键性质包括：对称性 P^T = P，幂等性 P^2 = P，以及 \text{rank}(P) = p+1。这些性质保证了 \mathrm{I}_n - P 也是对称幂等的，且 \text{rank}(\mathrm{I}_n - P) = n - p - 1，这解释了自由度的来源。

单个残差 \hat\varepsilon_i 的方差为：

\text{Var}(\hat\varepsilon_i) = \sigma^2(1 - P_{ii})

其中 P_{ii} 是投影矩阵的第 i 个对角元素，称为杠杆值（leverage）。杠杆值衡量了第 i 个观测点在设计空间中的"极端"程度。因此，标准化残差：

\tilde{\varepsilon}_i = \frac{\hat{\varepsilon}*i}{\sigma\sqrt{1-P*{ii}}}

服从标准正态分布。用 \hat\sigma 替代未知的 \sigma 得到学生化残差，它近似服从t分布。

残差诊断图

两个关键的诊断图帮助评估模型假设：

# Residual against fitted values
axes[0].scatter(OLSModel.fittedvalues, OLSModel.resid)
axes[0].axhline(y=0, color='grey', linestyle='dashed')

# QQ-plot of the residual
resStu = OLSModel.resid_pearson / np.sqrt(1 - OLSModel.get_influence().hat_matrix_diag)
stats.probplot(resStu, dist="norm", plot=axes[1])

残差vs拟合值图用于检测：异方差性（漏斗形模式），非线性关系（系统性模式），以及异常值。理想情况下，残差应随机分布在零线两侧，无明显模式。

Q-Q图（Quantile-Quantile plot）比较残差的经验分位数与理论正态分位数。如果残差服从正态分布，点应该落在45度直线上。常见的偏离模式包括：S形曲线（尾部过轻），反S形（尾部过重），以及离群点。

OLSModel.get_influence().hat_matrix_diag 返回杠杆值，用于计算学生化残差。高杠杆值的点对回归线有较大影响，需要特别关注。

练习3：部分系数估计（理论）

这个理论问题探讨了分组回归的重要技术。将回归变量分为两组 \underline{X}*1 和 \underline{X}\*2，对应系数为 \beta\*{(1)} 和 \beta*{(2)}：

\hat\beta_{(1)}, \hat\beta_{(2)} = \arg\min_{\beta_{(1)}, \beta_{(2)}} |\underline{y} - \underline{X}*1\beta*{(1)} - \underline{X}*2\beta*{(2)}|^2

如果只关心 \beta_{(2)} 而不想显式估计 \beta_{(1)}，可以使用Frisch-Waugh-Lovell定理。该定理指出，\hat\beta_{(2)} 可以通过以下步骤获得：

  1. 将 \underline{y} 对 \underline{X}_1 回归，得到残差 \underline{e}_y = (\mathrm{I} - P_1)\underline{y}
  2. 将 \underline{X}_2 的每一列对 \underline{X}*1 回归，得到残差矩阵 \underline{E}*{X_2} = (\mathrm{I} - P_1)\underline{X}_2
  3. 将 \underline{e}*y 对 \underline{E}*{X_2} 回归，得到的系数就是 \hat\beta_{(2)}

这个结果的直观解释是：\hat\beta_{(2)} 衡量的是 \underline{X}_2 中不能被 \underline{X}_1 解释的部分对 \underline{y} 中不能被 \underline{X}_1 解释的部分的影响。这种部分回归（partialing out）技术在控制混淆变量、理解条件关系以及计算效率方面都有重要应用。

从计算角度看，当 \underline{X}_1 包含大量控制变量时，这种方法可以避免求解高维线性系统。从因果推断角度看，它提供了在控制其他变量后评估特定变量效应的方法。

练习4：真实数据应用

California房价数据集介绍

California housing数据集是机器学习领域的经典数据集，包含了加州各地区的房价中位数及相关的社会经济和地理变量。这个数据集的特点包括：地理坐标（经度、纬度），人口统计（人口、家庭数），房屋特征（房间数、卧室数、房龄），以及经济指标（收入中位数、房价中位数）。

数据加载和初步探索：

import pandas as pd
df = pd.read_csv('california_housing_train.csv', sep=',')

# 查看数据基本信息
print(df.describe())  # 描述性统计
print(df.info())      # 数据类型和缺失值

# 散点图矩阵展示变量间关系
pd.plotting.scatter_matrix(df, alpha=0.8, figsize=(10, 10), diagonal='hist')

describe() 提供了每个变量的均值、标准差、分位数等统计信息，帮助识别潜在的异常值和了解数据分布。scatter_matrix 创建了所有变量两两之间的散点图，对角线显示每个变量的直方图，这是探索性数据分析的有力工具。

从散点图矩阵可以观察到几个重要模式：经纬度显示了加州的地理轮廓，某些变量对（如总房间数与家庭数）存在强相关性，房价分布呈现右偏，需要考虑对数变换，部分变量关系可能是非线性的。

多重共线性分析

多重共线性是指解释变量之间存在高度相关性，这会导致参数估计不稳定、标准误膨胀、以及解释困难。实验通过多种方法诊断多重共线性：

from statsmodels.stats.outliers_influence import variance_inflation_factor

# 计算方差膨胀因子（VIF）
vif_data = pd.DataFrame()
vif_data["Feature"] = X.columns
vif_data["VIF"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
vif_data = vif_data.sort_values('VIF', ascending=False)

VIF的计算公式为 \text{VIF}_j = \frac{1}{1-R_j^2}，其中 R_j^2 是将第 j 个变量对其他所有变量回归的决定系数。VIF的解释：

  • VIF = 1：无多重共线性
  • VIF < 5：轻微多重共线性，通常可接受
  • VIF 5-10：中度多重共线性，需要注意
  • VIF > 10：严重多重共线性，需要处理

相关矩阵热力图的绘制：

correlation_matrix = X.corr()
sns.heatmap(correlation_matrix, annot=True, fmt='.2f', cmap='coolwarm', 
            center=0, square=True, ax=axes[0,0])

cmap='coolwarm' 使用冷暖色调表示负相关和正相关，center=0 将颜色映射的中心设为0，使得无相关性呈现中性色。annot=True 在每个格子中显示相关系数的数值。

从分析结果可以看到，total_rooms、total_bedrooms、population 和 households 这四个变量存在严重的多重共线性（VIF > 100）。这是合理的，因为它们都是衡量区域规模的指标。经度和纬度的VIF也很高，但这可能反映了地理位置的非线性效应而非真正的共线性问题。

模型拟合与诊断

使用statsmodels进行回归分析，因为它提供了详细的统计输出：

# 添加截距并拟合模型
X_with_intercept = sm.add_constant(X)
housing_model = sm.OLS(y, X_with_intercept).fit()
print(housing_model.summary())

模型摘要包含了丰富的信息：

  • 系数估计与显著性：每个变量的系数、标准误、t统计量和p值
  • 模型拟合度：R²和调整R²衡量模型解释力
  • F统计量：检验所有系数（除截距外）是否同时为0
  • 信息准则：AIC、BIC用于模型比较
  • 诊断统计：Durbin-Watson检验序列相关性，Jarque-Bera检验残差正态性

残差诊断与模型改进

原始模型的残差诊断揭示了严重的问题：

# Residuals vs Fitted图显示异方差
axes[0].scatter(housing_model.fittedvalues, housing_model.resid)

# Q-Q图显示偏离正态分布
stats.probplot(resStu, dist="norm", plot=axes[1])

残差图呈现明显的喇叭形，表明异方差性：房价越高，预测误差的变异越大。这违反了同方差假设，会导致标准误估计偏误和假设检验失效。Q-Q图显示残差分布有厚尾，偏离正态分布。

对数变换是处理这些问题的常用方法：

# 对目标变量进行对数变换
y_log = np.log(y)  # 房价都是正值，直接取对数

# 重新拟合模型
housing_model_log = sm.OLS(y_log, X_reduced_with_const).fit()

对数变换的好处包括：

  1. 稳定方差：将乘性误差转换为加性误差
  2. 改善正态性：右偏分布经对数变换后更接近正态
  3. 经济学解释：系数可解释为弹性或百分比变化
  4. 缓解异常值影响：压缩大值的影响

变换后的模型诊断图显示了显著改善：残差的异方差性减弱，Q-Q图更接近直线，表明残差更接近正态分布。

模型解释与实践考虑

在对数模型中，系数的解释发生了变化。对于连续变量，系数表示自变量增加一个单位时，因变量的百分比变化（近似值）。例如，如果收入中位数的系数是0.5，意味着收入中位数增加1个单位，房价中位数增加约50%。

实践中还需要考虑的问题包括：

  1. 空间相关性：相邻地区的房价可能相关，违反独立性假设。可以考虑空间回归模型或聚类标准误。
  2. 非线性关系：某些变量（如房龄）可能与房价存在非线性关系。可以添加多项式项或使用样条函数。
  3. 交互效应：位置与其他特征的交互可能很重要。海边的房间数价值可能高于内陆。
  4. 缺失数据：实际应用中需要处理缺失值，可以使用删除、插值或多重插补等方法。
  5. 预测vs因果推断：这个模型适合预测，但不能直接用于因果推断。例如，增加房间数的系数不能解释为因果效应，因为可能存在未观测的混淆因素。
统计学习（整合版，第 1–4 章）
统计学习实验二：正则化回归与模型选择

评论