跳转至

牛顿法和拟牛顿 BFGS 法实现 Logistic 回归

推导二元 Logistic 回归的 Hessian 矩阵,利用牛顿法和拟牛顿 BFGS 法求回归系数的极大似然估计。所得模型在训练样本的预测准确度为 78%。

unnamed-chunk-7-1

问题描述:

problem-statement

数据文件:

data

推导似然函数的 Hessian 矩阵

似然函数的梯度

gradient

似然函数的 Hessian 矩阵

hessian

牛顿法算法步骤

newton-method

回溯线搜索的牛顿法

newton-method-backtracking

拟牛顿 BFGS 法算法步骤

BFGS-1

BFGS-2

拟牛顿 BFGS 法中\(B\)矩阵的迭代公式

BFGS-3

BFGS-4

使用 Python 将牛顿法和拟牛顿 BFGS 法应用于求解\(\beta\)的极大似然估计

Python
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
# 支持中文
matplotlib.rcParams['font.sans-serif'] = ['SimSun']
matplotlib.rcParams['axes.unicode_minus'] = False

导入数据

Python
df = pd.read_csv('data.csv')
X = df.iloc[:,2:].to_numpy()
Y = df.iloc[:,1:2].to_numpy()

求对数似然函数的相反数值,即\(J(\beta)\)

Python
def f(X, Y, beta):
    sum = 0
    for i in range(Y.shape[0]):
        product = np.exp(beta@X[i])
        if Y[i][0] == 1:
            sum = sum + np.log(product/(1+product))
        else:
            sum = sum + np.log(1/(1+product))
    return -sum

求梯度向量

Python
def gradient(X, Y, beta):
    # 将梯度向量初始化为零向量
    gradient_vector = np.zeros(X.shape[1])
    # 将 i 在 1 到 n 上遍历并加和,求梯度的值
    for i in range(Y.shape[0]):
        product = np.exp(beta@X[i])
        if Y[i][0] == 1:
            gradient_vector = gradient_vector - X[i] / (1 + product)
        else:
            gradient_vector = gradient_vector + product * X[i] / (1 + product)
    return gradient_vector

求 Hessian 矩阵。此题中的 Hessian 矩阵只与\(X\)\(\beta\)有关,与\(Y\)无关

Python
def hessian(X, beta):
    # 将 S 矩阵初始化。注意要初始化为 0.0,如果初始化为 0,则不能储存小数
    S = np.diag(np.full(Y.shape[0], 0.0))
    for i in range(Y.shape[0]):
        product = np.exp(beta@X[i])
        u = product / (1 + product)
        S[i, i] = u*(1 - u)
    return X.T@S@X

应用牛顿法,结合回溯线搜索求\(\beta\)的极大似然估计

值得注意的是,最初我没有用回溯线搜索控制步长,导致\(\beta\)的模长很大,代入到梯度公式中,由于要计算指数,而\(\beta\)出现在指数中,从而指数的结果过于巨大甚至溢出。这说明不控制步长的话结果可能不收敛,反而在朝目标值的反方向迭代。

由此可见,回溯线搜索等控制步长的方法具有很大的实际意义。

Python
# 设置 beta 的初始点
beta = np.ones(X.shape[1])
# 计算当前梯度的模
norm = np.linalg.norm(gradient(X, Y, beta))
norm_list = [norm]
# 当梯度的模长大于精度要求时,需要一直迭代
while norm > 1e-5:
    # 求 Hessian 矩阵
    hessian_matrix = hessian(X, beta)
    # 求 Hessian 矩阵的逆
    hessian_matrix_inverse = np.linalg.inv(hessian_matrix)
    # 求 beta 的迭代方向和长度
    v = - hessian_matrix_inverse@gradient(X, Y, beta)
    # 设置步长倍数的初始值
    t = 1
    # 回溯线性搜索,检查目前的步长是否满足要求,若步长太大,则需要缩小步长为原来的 0.5
    while f(X, Y, beta + t*v) > f(X, Y, beta) + 0.4 * t * gradient(X, Y, beta) @ v:
        t = 0.5 * t
    # 迭代更新 beta
    beta = beta + t*v
    # 求当前梯度的模长
    norm = np.linalg.norm(gradient(X, Y, beta))
    # 记录当前梯度的模长
    norm_list.append(norm)
print('beta 的最大似然估计是:', beta)
Text Only
## beta的最大似然估计是: [ 1.27459903  1.43703603  0.65622366 -1.22337458 -2.34681303]

绘制迭代次数与梯度的模

Python
fig = plt.figure(figsize=(5, 4), dpi=200)
plt.plot(norm_list, 'r-o', markerfacecolor='none')# 空心圆用 markerfacecolor='none',
plt.xlabel('迭代次数')
plt.ylabel('梯度的模')

unnamed-chunk-7-1

应用拟牛顿法,结合回溯线搜索求\(\beta\)的极大似然估计

Python
# 设置 beta 的初始点
beta = np.ones(X.shape[1])
# 设置 B 的初始值,为一个对称正定矩阵
B = np.diag(np.full(X.shape[1], 10))
# 当梯度的模长大于精度要求时,需要一直迭代
# 求 g_k
g_k = gradient(X, Y, beta)
# 计算当前梯度的模
norm = np.linalg.norm(g_k)
norm_list = [norm]
while norm > 1e-5:
    # 求迭代的方向和步长
    p = -np.linalg.inv(B)@g_k
    # 设置步长倍数的初始值
    t = 1
    # 回溯线性搜索,检查目前的步长是否满足要求,若步长太大,则需要缩小步长为原来的 0.5
    while f(X, Y, beta + t*p) > f(X, Y, beta) + 0.4 * t * g_k @ p:
        t = 0.5 * t
    # 迭代更新 B
    # 求 g_{k+1}
    g_k_plus_1 = gradient(X, Y, beta + t*p)
    # 求 y_k
    y_k = g_k_plus_1 - g_k
    # 求 sigma_k
    sigma_k = t*p
    # 迭代更新 B
    B = B + (np.outer(y_k,y_k))/(y_k@sigma_k) - (np.outer(B.dot(sigma_k),(sigma_k.T)).dot(B.T))/(sigma_k.T@B@sigma_k)
    # 迭代更新 beta
    beta = beta + t*p
    # 求当前梯度
    g_k = gradient(X, Y, beta)
    # 求当前梯度的模长
    norm = np.linalg.norm(g_k)
    # 记录当前梯度的模长
    norm_list.append(norm)
print('beta 的最大似然估计是:', beta)
Text Only
## beta的最大似然估计是: [ 1.27459939  1.43703625  0.65622384 -1.22337487 -2.34681354]

绘制迭代次数与梯度的模

Python
fig = plt.figure(figsize=(5, 4), dpi=200)
plt.plot(norm_list, 'r-o', markerfacecolor='none')# 空心圆用 markerfacecolor='none',
plt.xlabel('迭代次数', usetex=False) # 中文不用 tex 渲染,否则会报错
plt.ylabel('梯度的模', usetex=False)

unnamed-chunk-9-3

利用训练出的\(\beta\)值,对\(Y\)进行预测

Python
df['p_predict'] = np.nan
df['Y_predict'] = np.nan
df['correct'] = np.nan
df['p_predict'] = df.apply(lambda line: 1/(1+np.exp(-line.iloc[2:7]@beta)), axis=1)
# 设定判断阈值
hurdle = 0.5
df['Y_predict'] = df.apply(lambda line: 1 if line['p_predict']>hurdle else 0, axis=1)
df['correct'] = df.apply(lambda line: 1 if line['Y_predict']==line['Y'] else 0, axis=1)

计算预测准确率

Python
df
Text Only
##      Unnamed: 0  Y        X1        X2  ...        X5  p_predict  Y_predict  correct
## 0             1  0  0.227077 -0.492374  ...  0.337587   0.074140          0        1
## 1             2  1 -1.439121 -0.898333  ... -1.729708   0.840268          1        1
## 2             3  0 -0.252111 -0.274725  ... -0.705981   0.654752          1        0
## 3             4  1 -1.736983 -1.197634  ... -1.266615   0.471904          0        0
## 4             5  1 -1.081011 -0.283816  ... -0.937780   0.639866          1        1
## ..          ... ..       ...       ...  ...       ...        ...        ...      ...
## 195         196  1 -0.555757 -0.483450  ... -1.764113   0.993497          1        1
## 196         197  0 -0.581243 -0.590711  ...  0.485322   0.057910          0        1
## 197         198  0 -1.752883 -0.934694  ...  0.876816   0.004954          0        1
## 198         199  1  0.397383  0.231328  ...  0.902995   0.121900          0        0
## 199         200  1  1.127998  1.487247  ...  0.228712   0.963261          1        1
## 
## [200 rows x 10 columns]
Python
print('分类预测准确率为:{0:.2%}'.format(df['correct'].mean()))
Text Only
## 分类预测准确率为:78.50%

评论