NeuromatchAcademy课程内容：线性回归与最大似然估计

教程2：基于最大似然估计的线性回归

引言

在统计学和机器学习中，线性回归是最基础且广泛应用的模型之一。本教程将介绍如何使用最大似然估计(Maximum Likelihood Estimation, MLE)方法来拟合线性回归模型，这是对最小二乘法的一种概率视角的扩展。

1. 最大似然估计基础

1.1 高斯噪声模型

在传统线性回归中，我们假设数据由线性关系加上噪声组成：

y = θx + ε

其中ε是噪声项。当我们将噪声建模为高斯分布时，可以更全面地理解数据生成过程。高斯分布（正态分布）的概率密度函数为：

N(x; μ, σ²) = (1/√(2πσ²)) * exp(-(x-μ)²/(2σ²))

在简单线性回归中，我们通常假设噪声ε服从均值为0，方差为1的标准正态分布：ε ∼ N(0, 1)。

1.2 概率模型

基于高斯噪声假设，我们可以将响应变量y也视为随机变量：

y ∼ N(θx, 1)

这意味着给定x和参数θ，y的概率密度函数为：

p(y|x,θ) = (1/√(2π)) * exp(-(y-θx)²/2)

1.3 似然函数

似然函数L(θ|x,y)衡量了在给定数据(x,y)下，参数θ的可能性。对于单个数据点，似然函数等于概率密度函数：

L(θ|x,y) = p(y|x,θ)

2. 实现与可视化

2.1 生成模拟数据

我们首先生成一些模拟数据来演示概念：

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

np.random.seed(121)
theta = 1.2  # 真实参数值
n_samples = 30
x = 10 * np.random.rand(n_samples)  # 均匀分布在[0,10)的输入
noise = np.random.randn(n_samples)  # 标准正态噪声
y = theta * x + noise  # 生成输出数据

2.2 可视化概率密度

我们可以可视化p(y|x,θ=1.2)的概率密度：

def plot_density_image(x, y, theta, sigma=1, ax=None):
    """绘制给定x,θ,σ下y的概率密度图像"""
    if ax is None:
        fig, ax = plt.subplots()
    
    xmin, xmax = np.floor(np.min(x)), np.ceil(np.max(x))
    ymin, ymax = np.floor(np.min(y)), np.ceil(np.max(y))
    xx = np.linspace(xmin, xmax, 50)
    yy = np.linspace(ymin, ymax, 50)
    
    surface = np.zeros((len(yy), len(xx)))
    for i, x_i in enumerate(xx):
        surface[:, i] = stats.norm(theta * x_i, sigma).pdf(yy)
    
    ax.set(xlabel='x', ylabel='y')
    return ax.imshow(surface, origin='lower', aspect='auto', 
                    vmin=0, vmax=None, cmap='Wistia',
                    extent=[xmin, xmax, ymin, ymax])

# 绘制密度图像
fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(10, 5))
im = plot_density_image(x, y, 1.2, ax=ax1)
plt.colorbar(im, ax=ax1)
ax1.axvline(8, color='k')
ax1.set(title=r'p(y | x, $θ$=1.2)')

# 绘制x=8时的概率密度函数
ylim = ax1.get_ylim()
yy = np.linspace(ylim[0], ylim[1], 50)
ax2.plot(yy, stats.norm(theta * 8, 1).pdf(yy), color='orange', linewidth=2)
ax2.set(title=r'p(y|x=8, $θ$=1.2)', xlabel='y', ylabel='probability density')
plt.show()

2.3 实现似然函数

def likelihood(theta, x, y):
    """计算给定θ下数据(x,y)的似然值"""
    sigma = 1  # 假设噪声标准差为1
    return stats.norm(theta * x, sigma).pdf(y)

3. 最大似然估计

最大似然估计的核心思想是找到使似然函数最大的参数θ。在实际应用中，我们通常最大化对数似然函数，因为：

1. 对数函数是单调递增的，最大化对数似然等价于最大化似然
2. 对数转换将乘积变为求和，计算更稳定

3.1 对数似然函数

对于整个数据集的对数似然函数：

log L(θ) = Σ log p(y_i|x_i,θ) = -n/2 log(2π) - 1/2 Σ(y_i - θx_i)²

可以看到，最大化对数似然等价于最小化Σ(y_i - θx_i)²，这与最小二乘法一致。

3.2 寻找最优θ

我们可以通过以下步骤找到最大似然估计：

1. 定义似然函数
2. 在参数空间搜索使似然最大的θ值
3. 可视化似然函数随θ的变化

# 计算不同θ值的似然
theta_range = np.linspace(0, 2, 100)
likelihoods = [np.prod(likelihood(t, x, y)) for t in theta_range]

# 找到最大似然对应的θ
mle_theta = theta_range[np.argmax(likelihoods)]
print(f"最大似然估计的θ值: {mle_theta:.2f}")

# 可视化
plt.plot(theta_range, likelihoods)
plt.axvline(mle_theta, color='r', linestyle='--', label='MLE')
plt.xlabel('θ')
plt.ylabel('Likelihood')
plt.title('Likelihood Function')
plt.legend()
plt.show()

4. 讨论与扩展

1. 高斯噪声假设：当噪声不服从高斯分布时，MLE可能不是最优估计
2. 方差估计：本例假设噪声方差σ²=1，实际中可以将其作为参数一起估计
3. 多参数模型：对于多元线性回归，MLE方法可以自然扩展到多参数情况

5. 总结

本教程介绍了：

• 将线性回归建模为概率模型
• 高斯噪声假设下的似然函数
• 最大似然估计方法及其实现
• MLE与最小二乘法的等价性

最大似然估计为线性回归提供了概率论基础，使我们能够量化参数估计的不确定性，并为后续的模型比较和选择奠定了基础。