简介

我们将在此章节用 python 自己实现一遍以下几种模型：

• 线性回归（linear regression）
• 对数线性回归（log-linear regression）
• 对数几率回归（logistic regression，或译作逻辑回归）
• 线性判别分析（Linear Discriminant Analysis，简称 LDA）

在此之前我们需要一点点 python 的编程基础，会用以下的库：

• numpy：用于数值计算，包含各种矩阵计算的函数。
• matplotlib：用于数据可视化，将数据绘制成统计图。我们需要使用这个库来画图。
• pandas：非必需，用于数据分析，便于处理数据。不过我们用到的数据比较简单，不需要这个库也可以处理样本数据。

线性回归（linear regression）

在中学的时候我们学过线性函数（一次函数），即

我们可以从上式推出对于多元线性函数的一般式，即

其中，x = （x1；x2；...；xn）。若令 k = （k1；k2；...；kn），则

再令 θ = （b；k），x^ = （1；x），则得到 f(x) 的一般式

其中 θ 和 x^ 均为 n + 1 维向量。对于每一个 n 维特征的样本，都可以用 n + 1维向量 x^ 表示，再用 f(x) 的一般式计算出其预测值。此时，我们只需要知道 θ 的值就能写出线性回归了。我们知道，线性模型是从样本中训练出来的，所以我们需要用样本集 D 来训练出 θ 的值。

令 m 表示样本集 D 的样本数量，d 表示每个样本的特征数量，则我们可以把样本集 D 表示为一个 m x (d + 1) 大小的矩阵 X，即

令向量 y 表示 D 的标签（label），即

通过第二章的学习，我们知道均方误差是回归任务中最常用的性能度量。而最小化均方误差即能使线性模型 f(x) 获得最好的性能，所以 θ 应是令 f(x) 的均方误差最小化的值，即

令

对 θ^ 求导得

令上式为零可得 θ^ 的最优解的闭式（closed-form）解，当 XTX为满秩矩阵（full-rank matrix）或正定矩阵（positive definite matrix）时，令上式为零可得

即得到 θ 的最优解。此方法被称为最小二乘法（least square method），即基于均方误差最小化来进行模型求解的方法。

最终得到线性回归模型：

线性回归（Linear Regression）

以下为线性回归的 python 实现：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

print('正在读取数据...')

data = pd.read_csv('data_regression.csv')

TRAIN_NUM = len(data)
TRAIN_DIM = data.ndim - 1
train_data = data.values[0:TRAIN_NUM, 0:TRAIN_DIM]
train_label = data.values[0:TRAIN_NUM, TRAIN_DIM]

# 处理训练集数据
# TRAIN_NUM   : m
# TRAIN_DIM   : n
# train_data  : m x n
# train_label : 1 x m

print('数据读取完成，总共 ', TRAIN_NUM, ' 行数据')

def linear(data_x, data_y):
    # data_x : train_data
    # data_y : train_label
    # x      : [1] + train_data, m x n
    # y      : train_label.T, m x 1
    m = len(data_y)
    n = np.size(data_x, 1) + 1
    x = np.mat([[1] + list(data_x[i]) for i in range(m)])
    y = np.mat(data_y).T
    theta = (x.T * x).I * x.T * y
    return theta

print('正在计算 Theta 值...')

theta = linear(train_data, train_label)

print('计算 Theta 值成功')
print('Theta = ', theta)

n_grids = 100
x_axis = np.linspace(1, 100, n_grids)
y_axis = theta.T * np.mat([[1] + [x_axis[i]] for i in range(n_grids)]).T
plt.figure('Linear Regression')
plt.plot(x_axis, y_axis.tolist()[0], 'k')
plt.scatter(train_data[0:TRAIN_NUM, 0:TRAIN_DIM], train_label)
plt.show()

数据文件 data_regression.csv

data,label
58,187.8
2.7,26.72
99.4,322.84
10.1,68.36
22.7,118.72
39.9,101.64
97.4,399.64
51,226.6
14.6,25.56
25.7,65.52
27,84.2
26.5,144.4
11.8,21.48
11.1,34.96
69,207.4
14,19.4
69.4,283.84
35.6,93.16
14.8,19.28
51.7,179.12

另外，将代码中计算的 theta 一行

theta = (x.T * x).I * x.T * y

改为

theta = (x.T * x).I * x.T * np.log(y)

再将计算图表 Y 轴的 y_axis 一行

y_axis = theta.T * np.mat([[1] + [x_axis[i]] for i in range(n_grids)]).T

改为

y_axis = np.exp(theta.T * np.mat([[1] + [x_axis[i]] for i in range(n_grids)]).T)

则将线性回归变为了对数线性回归（log-linear regression），即

对数线性回归（Log-Linear Regression）

更一般地，考虑单调可微函数 g（·），令

我们可以得到广义线性模型（generalized linear model），其中函数 g（·） 称为联系函数（link function）。

显然，上面的对数线性回归就是在 g（·） = ln（·）时的特例。

对数几率回归（logistic regression，或逻辑回归）

对数几率回归与线性回归相似，主要区别在于对数几率回归使用了 Sigmoid 函数，即形似 S 的函数，如

将线性回归的一般式代入 z，可得

其中，y ∈ （0，1）。

对于二分类任务，可以令函数 y 代表样本被分类为正例的概率，令函数 1 - y 代表样本被分类为反例的概率，则

其中 p（y = i | x）表示数据被分类为 i 的概率。二分类任务中，1 表示正例，0 表示反例。此时，我们就可以通过极大似然法（maximum likelihood method）来估计 θ 值，则有

要使 θ 得到最优解，我们就要最大化 ℓ(θ)，即最小化 -ℓ(θ)，则有

要使 θ 达到最优解，可用梯度下降法（gradient descent method）、牛顿法（Newton method）等。我们用牛顿法来求其最优解，其第 t + 1 轮的更新公式为

其中，关于 θ 的一阶导数、二阶导数分别为

注意在矩阵运算中的 p（y = 1 | x）是 m x 1 的矩阵，而在级数中的 p（y = i | x）是数（即 1 x 1 的矩阵）

最终得到对数几率回归模型：

对数几率回归（Logistic Regression）

以下为对数几率回归的 python 实现：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

print('正在读取数据...')

data = pd.read_csv('data_classification.csv')

TRAIN_NUM = len(data)
TRAIN_DIM = data.ndim - 1
train_data = data.values[0:TRAIN_NUM, 0:TRAIN_DIM]
train_label = data.values[0:TRAIN_NUM, TRAIN_DIM]

# 处理训练集数据
# TRAIN_NUM   : m
# TRAIN_DIM   : n
# train_data  : m x n
# train_label : 1 x m

print('数据读取完成，总共', TRAIN_NUM, '行数据')

def p1(x, theta):
    # x     : m x n
    # theta : n x 1
    return 1 - 1 / (1 + np.exp(x * theta))

def p1_single(x, theta):
    # x     : 1 x n
    # theta : n x 1
    return 1 - 1 / (1 + np.exp(theta.T * x.T))

def newtonsMethod(x, y, theta):
    # x     : m x n
    # y     : m x 1
    # theta : n x 1
    m = len(y)
    n = np.size(x, 1)
    fir_d = (1 / m) * x.T * (p1(x, theta) - y)
    sec_d = 0
    for i in range(m):
        sec_d += x[i].T * x[i] * float(p1_single(x[i], theta) * (1 - p1_single(x[i], theta)))
    sec_d /= m
    return theta - sec_d.I * fir_d

def costFunction(x, y, theta):
    # x     : m x n
    # y     : m x 1
    # theta : n x 1
    m = len(y)
    return -1 / m * (y.T * np.log(p1(x, theta)) + (1 - y).T * np.log(1 - p1(x, theta)))

def logistic(data_x, data_y, max_times):
    # data_x    : train_data
    # data_y    : train_label
    # max_times : int
    # x         : [1] + train_data, m x n
    # y         : train_label.T, m x 1
    m = len(data_y)
    n = np.size(data_x, 1) + 1
    x = np.mat([[1] + list(data_x[i]) for i in range(m)])
    y = np.mat(data_y).T
    theta = np.zeros((n, 1))
    times = 0
    J = costFunction(x, y, theta)
    while True:
        last_J = J
        theta = newtonsMethod(x, y, theta)
        J = costFunction(x, y, theta)
        times += 1
        if (last_J - J < 1e-6) or times >= max_times:
            break
    print('总共迭代', times, '次')
    return theta

print('正在计算 Theta 值...')

theta = logistic(train_data, train_label, 20)

print('计算 Theta 值成功')
print('Theta = ', theta)

n_grids = 100
x_axis = np.linspace(0, 100, n_grids)
y_axis = p1(np.mat([[1] + [x_axis[i]] for i in range(n_grids)]), theta)
plt.figure('Linear Regression')
plt.plot(x_axis, y_axis, 'k')
plt.scatter(train_data[0:TRAIN_NUM, 0:TRAIN_DIM], train_label)

y_label = (p1(np.mat([[1] + [train_data[0:TRAIN_NUM, 0:TRAIN_DIM][i]] for i in range(TRAIN_NUM)]), theta) > 0.5).T * 1
correct = np.sum((abs(y_label[0] - train_label) < 1e-3) * 1)
print('分类正确率为：', correct / TRAIN_NUM * 100, '%')
plt.scatter(train_data[0:TRAIN_NUM, 0:TRAIN_DIM], y_label)

plt.show()

数据文件 data_classification.csv

data,label
92.44,0
82,0
40.47,0
25.29,1
18.71,1
51.88,0
69.92,0
29.58,1
50.5,1
65.59,0
81.42,0
89.93,0
89.08,0
17.11,1
20.75,0
22.51,0
77.43,1
60.48,0
41.47,1
54.86,0

线性判别分析（Linear Discriminant Analysis，简称 LDA）

公式太多，暂时推不过来，有缘再加……

小结

这一章给出了几个基本的线性模型，并且介绍了下二分类学习及多分类学习的方法，优化算法方面着重讲解了最小二乘法和牛顿法，却没有细讲很是通用的梯度下降法。

最小二乘法是用矩阵计算快速给出 θ 的最优解，不过这种方法对于特征数量很大的样本集来说过于缓慢，几乎没法计算。

牛顿法相较于梯度下降法来说，下降的幅度更大，能用更少的迭代次数给出最优解。不过和最小二乘法相似的是，牛顿法对于特征数量大的样本集来说也过于缓慢，主要是由于在求二阶导数的时候用到了黑塞矩阵（Hessian Matrix）。对于此类问题，选择拟牛顿法（Quasi-Newton Methods）更加合适。