简介
我们将在此章节用 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），即基于均方误差最小化来进行模型求解的方法。
最终得到线性回归模型：
以下为线性回归的 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），即
更一般地，考虑单调可微函数 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 的矩阵）
最终得到对数几率回归模型：
以下为对数几率回归的 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）更加合适。
