6统计学习方法(第2版)

逻辑斯谛回归与最大熵模型

对数线性模型

逻辑斯谛回归最大熵模型改进的迭代尺度法

🎯 AI 时代的学习建议:不要只读理论。理解这一章的设计动机(为什么需要这个概念?)和边界条件(什么情况下会失效?)。当 AI 生成的代码涉及这些概念时,你能判断它是否用对了场景。AI 能帮你总结内容,但无法替代你亲手阅读、思考、建立直觉的过程。

第6章:逻辑斯谛回归与最大熵模型(Logistic Regression and Maximum Entropy Models)

一、导读

1.1 本章学习目标

第6章介绍逻辑斯谛回归——虽然名字里有"回归",但它实际上是一种分类算法。逻辑斯谛回归是最经典的线性分类方法之一,在工业界有着广泛的应用。同时,本章还介绍最大熵模型——一种基于信息论原理的通用分类框架。

本章的核心学习目标包括:

  • 理解逻辑斯谛回归模型:sigmoid函数如何将线性输出映射为概率
  • 掌握对数几率的概念:为什么 log(P/(1-P)) = w·x 是一个好的模型?
  • 理解学习策略:为什么用对数似然损失而不是平方损失?
  • 掌握优化算法:梯度下降法、牛顿法、拟牛顿法(BFGS)
  • 理解最大熵原理:为什么在满足约束的条件下要使熵最大?
  • 实现逻辑斯谛回归:从零开始编码,对比不同优化算法

1.2 为什么本章如此重要

逻辑斯谛回归是工业界最常用的分类算法之一。它的优势在于:简单、高效、可解释、概率输出。在广告推荐、信用评分、疾病预测等领域,逻辑斯谛回归往往是基线模型。

逻辑斯谛回归也是理解神经网络的起点。单个逻辑斯谛回归单元就是神经网络中最基本的"神经元"。理解逻辑斯谛回归,就是理解深度学习的基石。

1.3 核心问题引导

在阅读本章时,请带着以下问题思考:

逻辑斯谛回归的sigmoid函数有什么特殊性质?为什么选择它?

对数几率(log-odds)的直觉含义是什么?

为什么逻辑斯谛回归用对数似然损失而不是平方损失?

梯度下降法和牛顿法在优化逻辑斯谛回归时有什么区别?

最大熵模型和逻辑斯谛回归有什么关系?

逻辑斯谛回归如何处理多分类问题?


二、核心概念详解

6.1 逻辑斯谛回归模型

6.1.1 Sigmoid函数

逻辑斯谛回归的核心是sigmoid函数:

σ(z) = 1 / (1 + e^(-z))

sigmoid函数将任意实数映射到 (0, 1) 区间,可以解释为概率。

python
import numpy as np

def sigmoid(z):
    return 1 / (1 + np.exp(-np.clip(z, -500, 500)))

# sigmoid函数的性质
z = np.linspace(-10, 10, 100)
print(f"sigmoid(0) = {sigmoid(0):.4f}")
print(f"sigmoid(∞) → {sigmoid(100):.4f}")
print(f"sigmoid(-∞) → {sigmoid(-100):.4f}")
print(f"sigmoid'(0) = {sigmoid(0) * (1 - sigmoid(0)):.4f}")
print("\nsigmoid的关键性质:")
print("1. 输出范围 (0, 1),可以解释为概率")
print("2. 单调递增")
print("3. 导数 σ'(z) = σ(z)(1-σ(z)),计算方便")
print("4. z=0时输出0.5,是决策边界")

6.1.2 模型定义

二分类逻辑斯谛回归模型:

P(Y=1|x) = 1 / (1 + e^(-(w·x + b)))
P(Y=0|x) = 1 - P(Y=1|x) = 1 / (1 + e^(w·x + b))

决策规则:如果 P(Y=1|x) ≥ 0.5,则预测为正类。

6.1.3 对数几率

逻辑斯谛回归的一个优美性质是:

log(P(Y=1|x) / P(Y=0|x)) = w·x + b

即对数几率(log-odds)是输入的线性函数。这意味着:

  • w·x 每增加1,几率(odds)变为原来的 e 倍
  • w_i 表示特征 x_i 每增加1个单位,对数几率的变化量
python
import numpy as np

def log_odds_interpretation():
    """对数几率的直观解释"""
    w = np.array([0.5, -0.3, 1.0])
    b = -0.5
    
    # 基准点
    x0 = np.array([0, 0, 0])
    z0 = w @ x0 + b
    p0 = sigmoid(z0)
    
    print(f"基准点 x={x0}: P(Y=1) = {p0:.4f}")
    print()
    
    # 每个特征增加1个单位
    for i in range(3):
        x1 = x0.copy()
        x1[i] = 1
        z1 = w @ x1 + b
        p1 = sigmoid(z1)
        odds_ratio = (p1 / (1 - p1)) / (p0 / (1 - p0))
        print(f"特征{i}增加1: P(Y=1) = {p1:.4f}, 几率比 = {odds_ratio:.4f} = e^{w[i]:.1f}")
    
    print(f"\nw[0]=0.5: 特征0增加1,对数几率增加0.5")
    print(f"w[1]=-0.3: 特征1增加1,对数几率减少0.3")
    print(f"w[2]=1.0: 特征2增加1,对数几率增加1.0(影响最大)")

log_odds_interpretation()

6.2 学习策略:对数似然损失

6.2.1 为什么不用平方损失

逻辑斯谛回归不使用平方损失,原因有两个:

非凸性:平方损失在逻辑斯谛回归中不是凸函数,存在局部最优

不匹配:平方损失假设误差是高斯分布的,但分类问题的输出是伯努利分布

6.2.2 对数似然损失

逻辑斯谛回归使用对数似然损失(Log Loss / Cross-Entropy Loss):

L(w, b) = -(1/N) * Σ_i [y_i * log(p_i) + (1-y_i) * log(1-p_i)]

其中 p_i = P(Y=1|x_i) = sigmoid(w·x_i + b)。

这个损失函数是凸函数,有唯一的全局最优解。

python
import numpy as np

def log_loss(y_true, y_pred):
    """对数似然损失"""
    eps = 1e-15
    y_pred = np.clip(y_pred, eps, 1 - eps)
    return -np.mean(y_true * np.log(y_pred) + (1 - y_true) * np.log(1 - y_pred))

def squared_loss(y_true, y_pred):
    """平方损失"""
    return np.mean((y_true - y_pred) ** 2)

# 对比两种损失函数的梯度行为
y_true = np.array([1, 0, 1, 0])
y_pred_wrong = np.array([0.1, 0.9, 0.1, 0.9])  # 错误预测
y_pred_right = np.array([0.9, 0.1, 0.9, 0.1])  # 正确预测

print(f"{'':>15} {'错误预测':>12} {'正确预测':>12}")
print("-" * 42)
print(f"{'对数损失':>15} {log_loss(y_true, y_pred_wrong):>12.4f} {log_loss(y_true, y_pred_right):>12.4f}")
print(f"{'平方损失':>15} {squared_loss(y_true, y_pred_wrong):>12.4f} {squared_loss(y_true, y_pred_right):>12.4f}")

6.3 学习算法

6.3.1 梯度下降法

逻辑斯谛回归的梯度:

∂L/∂w = (1/N) * X^T · (p - y)
∂L/∂b = (1/N) * Σ(p_i - y_i)

其中 p = sigmoid(X·w + b)。

python
import numpy as np

class LogisticRegressionGD:
    def __init__(self, lr=0.01, n_iter=1000):
        self.lr = lr
        self.n_iter = n_iter
        self.loss_history = []
    
    def fit(self, X, y):
        self.w = np.zeros(X.shape[1])
        self.b = 0.0
        
        for _ in range(self.n_iter):
            z = X @ self.w + self.b
            p = sigmoid(z)
            
            # 梯度
            dw = X.T @ (p - y) / len(y)
            db = np.mean(p - y)
            
            # 更新
            self.w -= self.lr * dw
            self.b -= self.lr * db
            
            # 记录损失
            loss = log_loss(y, p)
            self.loss_history.append(loss)
        
        return self
    
    def predict_proba(self, X):
        return sigmoid(X @ self.w + self.b)
    
    def predict(self, X):
        return (self.predict_proba(X) >= 0.5).astype(int)

# 测试
from sklearn.datasets import make_classification
X, y = make_classification(n_samples=500, n_features=10, random_state=42)
X = (X - X.mean(axis=0)) / X.std(axis=0)

model = LogisticRegressionGD(lr=0.1, n_iter=500)
model.fit(X, y)
print(f"梯度下降法最终损失: {model.loss_history[-1]:.6f}")
print(f"训练准确率: {np.mean(model.predict(X) == y):.4f}")

6.3.2 牛顿法

牛顿法利用二阶导数(Hessian矩阵)信息,收敛速度更快:

w_{t+1} = w_t - H^{-1} · g

其中 g 是梯度,H 是 Hessian 矩阵。

python
import numpy as np

class LogisticRegressionNewton:
    def __init__(self, n_iter=20):
        self.n_iter = n_iter
        self.loss_history = []
    
    def fit(self, X, y):
        n, p = X.shape
        self.w = np.zeros(p)
        self.b = 0.0
        
        # 合并偏置
        X_aug = np.column_stack([X, np.ones(n)])
        theta = np.zeros(p + 1)
        
        for _ in range(self.n_iter):
            z = X_aug @ theta
            pred = sigmoid(z)
            
            # 梯度
            grad = X_aug.T @ (pred - y) / n
            
            # Hessian矩阵
            W = np.diag(pred * (1 - pred))
            H = X_aug.T @ W @ X_aug / n
            
            # 牛顿更新
            theta -= np.linalg.solve(H + 1e-8 * np.eye(p + 1), grad)
            
            loss = log_loss(y, pred)
            self.loss_history.append(loss)
        
        self.w = theta[:-1]
        self.b = theta[-1]
        return self
    
    def predict_proba(self, X):
        return sigmoid(X @ self.w + self.b)
    
    def predict(self, X):
        return (self.predict_proba(X) >= 0.5).astype(int)

# 对比收敛速度
model_gd = LogisticRegressionGD(lr=0.1, n_iter=500)
model_gd.fit(X, y)

model_newton = LogisticRegressionNewton(n_iter=20)
model_newton.fit(X, y)

print(f"梯度下降法: {len(model_gd.loss_history)}次迭代, 最终损失={model_gd.loss_history[-1]:.6f}")
print(f"牛顿法: {len(model_newton.loss_history)}次迭代, 最终损失={model_newton.loss_history[-1]:.6f}")
print(f"牛顿法收敛快得多!但每次迭代计算Hessian矩阵的代价更高")

6.3.3 拟牛顿法(BFGS)

BFGS是一种拟牛顿法,它通过近似Hessian矩阵的逆来避免直接计算Hessian:

python
import numpy as np
from scipy.optimize import minimize

class LogisticRegressionBFGS:
    def fit(self, X, y):
        n, p = X.shape
        
        def objective(theta):
            w, b = theta[:-1], theta[-1]
            z = X @ w + b
            p_pred = sigmoid(z)
            loss = -np.mean(y * np.log(p_pred + 1e-10) + (1 - y) * np.log(1 - p_pred + 1e-10))
            return loss
        
        def gradient(theta):
            w, b = theta[:-1], theta[-1]
            z = X @ w + b
            p_pred = sigmoid(z)
            grad_w = X.T @ (p_pred - y) / n
            grad_b = np.mean(p_pred - y)
            return np.concatenate([grad_w, [grad_b]])
        
        theta0 = np.zeros(p + 1)
        result = minimize(objective, theta0, jac=gradient, method='BFGS')
        
        self.w = result.x[:-1]
        self.b = result.x[-1]
        return self
    
    def predict_proba(self, X):
        return sigmoid(X @ self.w + self.b)
    
    def predict(self, X):
        return (self.predict_proba(X) >= 0.5).astype(int)

model_bfgs = LogisticRegressionBFGS()
model_bfgs.fit(X, y)
print(f"BFGS准确率: {np.mean(model_bfgs.predict(X) == y):.4f}")

6.4 最大熵模型

6.4.1 最大熵原理

最大熵原理:在满足所有已知约束的条件下,选择熵最大的概率分布。

直觉:如果我们只知道部分信息,那么最"公平"的假设是让未知的部分尽可能均匀分布。

6.4.2 最大熵模型

最大熵模型是最大熵原理在分类问题中的应用。给定特征函数 f_i(x, y),最大熵模型的形式为:

P(y|x) = (1/Z(x)) * exp(Σ_i w_i * f_i(x, y))

其中 Z(x) 是归一化因子。

当特征函数是线性函数时,最大熵模型等价于逻辑斯谛回归。

python
import numpy as np

def max_entropy_demo():
    """最大熵原理的直观演示"""
    # 假设一个骰子,我们知道它的平均值是3.5
    # 最大熵分布是什么?
    
    # 约束:Σ p_i = 1, Σ i*p_i = 3.5
    # 最大熵分布:均匀分布 p_i = 1/6
    
    # 验证:均匀分布的均值
    p_uniform = np.ones(6) / 6
    mean_uniform = np.sum(np.arange(1, 7) * p_uniform)
    entropy_uniform = -np.sum(p_uniform * np.log2(p_uniform))
    
    print(f"均匀分布:")
    print(f"  均值 = {mean_uniform:.2f}")
    print(f"  熵 = {entropy_uniform:.4f} bits")
    
    # 对比一个非均匀分布
    p_biased = np.array([0.1, 0.1, 0.1, 0.3, 0.2, 0.2])
    mean_biased = np.sum(np.arange(1, 7) * p_biased)
    entropy_biased = -np.sum(p_biased * np.log2(p_biased))
    
    print(f"\n非均匀分布:")
    print(f"  均值 = {mean_biased:.2f}")
    print(f"  熵 = {entropy_biased:.4f} bits")
    
    print(f"\n均匀分布的熵更大,这就是最大熵原理")

max_entropy_demo()

6.5 多分类扩展:Softmax回归

逻辑斯谛回归可以自然地扩展到多分类问题,使用softmax函数:

P(Y=k|x) = exp(w_k·x) / Σ_j exp(w_j·x)
python
import numpy as np

def softmax(z):
    """softmax函数"""
    z_max = np.max(z, axis=-1, keepdims=True)
    exp_z = np.exp(z - z_max)
    return exp_z / np.sum(exp_z, axis=-1, keepdims=True)

class SoftmaxRegression:
    def __init__(self, lr=0.01, n_iter=500):
        self.lr = lr
        self.n_iter = n_iter
    
    def fit(self, X, y):
        n, p = X.shape
        self.classes = np.unique(y)
        K = len(self.classes)
        
        self.W = np.zeros((p, K))
        self.b = np.zeros(K)
        
        # one-hot编码
        Y_onehot = np.zeros((n, K))
        for i, c in enumerate(self.classes):
            Y_onehot[y == c, i] = 1
        
        for _ in range(self.n_iter):
            z = X @ self.W + self.b
            P = softmax(z)
            
            # 梯度
            error = P - Y_onehot
            self.W -= self.lr * X.T @ error / n
            self.b -= self.lr * np.mean(error, axis=0)
        
        return self
    
    def predict_proba(self, X):
        return softmax(X @ self.W + self.b)
    
    def predict(self, X):
        P = self.predict_proba(X)
        return self.classes[np.argmax(P, axis=1)]

# 测试
from sklearn.datasets import load_iris
iris = load_iris()
X = (iris.data - iris.data.mean(axis=0)) / iris.data.std(axis=0)
y = iris.target

sm = SoftmaxRegression(lr=0.1, n_iter=1000)
sm.fit(X, y)
print(f"Softmax回归准确率: {np.mean(sm.predict(X) == y):.4f}")

三、实战案例

3.1 正则化逻辑斯谛回归

python
import numpy as np

class RegularizedLogisticRegression:
    def __init__(self, lr=0.01, n_iter=1000, reg_type='l2', C=1.0):
        self.lr = lr
        self.n_iter = n_iter
        self.reg_type = reg_type
        self.C = C
    
    def fit(self, X, y):
        n, p = X.shape
        self.w = np.zeros(p)
        self.b = 0.0
        
        for _ in range(self.n_iter):
            z = X @ self.w + self.b
            pred = sigmoid(z)
            
            # 梯度
            dw = X.T @ (pred - y) / n
            db = np.mean(pred - y)
            
            # 正则化
            if self.reg_type == 'l2':
                dw += self.w / self.C
            elif self.reg_type == 'l1':
                dw += np.sign(self.w) / self.C
            
            self.w -= self.lr * dw
            self.b -= self.lr * db
        
        return self
    
    def predict(self, X):
        return (sigmoid(X @ self.w + self.b) >= 0.5).astype(int)

# 对比不同正则化
np.random.seed(42)
X, y = make_classification(n_samples=200, n_features=20, n_informative=5, random_state=42)
X = (X - X.mean(axis=0)) / X.std(axis=0)

for reg in ['none', 'l1', 'l2']:
    if reg == 'none':
        model = LogisticRegressionGD(lr=0.1, n_iter=1000)
    else:
        model = RegularizedLogisticRegression(lr=0.1, n_iter=1000, reg_type=reg, C=0.1)
    model.fit(X, y)
    acc = np.mean(model.predict(X) == y)
    n_nonzero = np.sum(np.abs(model.w) > 0.01)
    print(f"{reg:>5}: 准确率={acc:.4f}, 非零系数={n_nonzero}")

3.2 决策边界可视化

python
import numpy as np

def visualize_decision_boundary():
    """可视化逻辑斯谛回归的决策边界"""
    np.random.seed(42)
    
    # 生成二分类数据
    X_pos = np.random.randn(50, 2) + np.array([1.5, 1.5])
    X_neg = np.random.randn(50, 2) + np.array([-1, -1])
    X = np.vstack([X_pos, X_neg])
    y = np.array([1] * 50 + [0] * 50)
    
    # 训练
    model = LogisticRegressionGD(lr=0.1, n_iter=500)
    model.fit(X, y)
    
    print(f"w = [{model.w[0]:.4f}, {model.w[1]:.4f}]")
    print(f"b = {model.b:.4f}")
    print(f"决策边界: {model.w[0]:.4f}*x1 + {model.w[1]:.4f}*x2 + {model.b:.4f} = 0")
    print(f"即 x2 = {-model.w[0]/model.w[1]:.4f}*x1 + {-model.b/model.w[1]:.4f}")
    print(f"准确率: {np.mean(model.predict(X) == y):.4f}")

visualize_decision_boundary()

四、本章小结

第6章介绍了逻辑斯谛回归与最大熵模型:

逻辑斯谛回归模型:通过sigmoid函数将线性输出映射为概率,对数几率是输入的线性函数。

学习策略:使用对数似然损失(交叉熵损失),它是凸函数,有唯一全局最优。

优化算法:梯度下降法简单但收敛慢;牛顿法利用二阶信息收敛快;BFGS拟牛顿法在实际中最常用。

最大熵模型:在满足约束条件下使熵最大,当特征为线性时等价于逻辑斯谛回归。

多分类扩展:通过softmax函数自然地扩展到多分类问题。

逻辑斯谛回归的核心启示是:简单的线性模型 + 非线性映射 = 强大的分类器。


五、思考与练习

思考题

sigmoid函数有什么特殊性质使得它特别适合分类问题?还有其他类似的函数吗?

为什么逻辑斯谛回归不用平方损失?从凸性和概率分布的角度解释。

牛顿法为什么比梯度下降法收敛快?它的缺点是什么?

最大熵原理的直觉是什么?为什么"最不确定"的分布反而最好?

逻辑斯谛回归在什么情况下会欠拟合?如何解决?

编程练习

从零实现逻辑斯谛回归(梯度下降法),在二分类数据上测试。

实现牛顿法优化逻辑斯谛回归,对比与梯度下降法的收敛速度。

实现Softmax回归,在鸢尾花数据集(三分类)上测试。


6/12逻辑斯谛
0%