机器学习 - 线性回归与逻辑回归（实践部分）

之前对线性回归和逻辑回归的理论部分做了较为详细的论述，下面通过一些例子再来巩固一下之前所学的内容。
需要说明的是，虽然我们在线性回归中都是直接通过公式推导求出w和b的精确值，但在实际运用中基本上都会采用梯度下降法作为首选，因为用代码表示公式会比较繁琐，而梯度下降法只需要不断对参数更新公式进行迭代即可，用代码表示非常简单，且迭代次数基本上不会对性能有太大的拖累，所以下面我们也将采用梯度下降的方法。

使用线性回归模拟数据分布

了解数据初始分布情况

现有一个数据集（提取码：yim9），描述了一系列点的分布情况。我们先将初始的分布情况绘制出来，便于对数据有个整体的了解：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

data = np.loadtxt('./data/linear_regression_data.txt', delimiter=',')
print(data[:10])  # 查看部分样本数据
data.shape        # 查看样本规模

# 将数据集转换成 y = w*x 的向量形式
x = np.c_[np.ones(data.shape[0]), data[:, 0]]  # 使用np.c_转换为列向量
y = np.c_[data[:, 1]]  # 1是给b留的位置

# 展示样本分布情况
plt.scatter(x[:, 1], y, color='green', marker='x', linewidth=1)
plt.xlabel('X')
plt.ylabel('Y')

可以看到，点在X轴左侧分布得比较密集，但是越往右，越呈现出按照直线分布的趋势，所以可以尝试使用线性回归进行模拟。

使用最小二乘法构造目标函数

我们需要求出目标函数，以便求出梯度，并记录目标函数在每轮迭代中的值

# 目标函数
def target_func(x, y, w=[[0],[0]]):
    # 这里我们设置w和b的初始值为0
    h = x.dot(w)                                   # h表示预测值
    m = y.size                                     # 通常会使用样本大小m平均整体误差的值
    target = 1.0/(2*m) * (np.sum(np.square(h-y)))  # 使用最小二乘法得到目标函数
    return target

计算梯度及参数更新公式

# 计算梯度
def grad(x, y, w):
    h = x.dot(w)
    m = y.size
    grad = (1.0/m)*(x.T.dot(h-y))  # 对目标函数求导即是梯度
    return grad

# 使用梯度循环更新参数
def grad_descent(x, y, nums_iter=1000, n=0.01, w=[[0],[0]]):
    # 默认迭代1000次，学习率为0.01
    target_his = np.zeros(nums_iter)    # 存放目标函数每轮迭代后的值
    for i in np.arange(nums_iter):
        # 循环更新w，并记录目标函数的值
        w = w - n*grad(x, y, w)
        target_his[i] = target_func(x, y, w)
    return w, target_his

w, target_his = grad_descent(x, y)

# 展示迭代效果
plt.plot(target_his)
plt.xlabel('Iter Times')
plt.ylabel('Loss Function')

可以看出，迭代1000次，性能逐渐趋于平稳。
最后与sklearn的结果相对比：

# 样本
plt.scatter(x[:, 1], y, color='red', marker='x', linewidth=1)
# 自主
xx = np.arange(5, 23)
yy = w[0] + xx*w[1]
plt.plot(xx, yy, label='LD(GR)')

# sklearn
regr = LinearRegression()
# reshape转换成列向量
regr.fit(x[:, 1].reshape(-1, 1), y.ravel())
yy2 = regr.intercept_ + regr.coef_ * xx
plt.plot(xx, yy2, label='LD(sklearn)')

plt.xlim(4, 24)
plt.xlabel('X')
plt.ylabel('Y')
plt.legend(loc=4)

可以看出效果基本一致。

使用逻辑回归分类

现有一个考试成绩数据集(提取码：7n04)，该数据集为某班级英语和数学两个科目的考试成绩，前两列分别表示英语和数学成绩，最后一列表示最终是否通过。

加载数据，查看数据样本及可视化

import pandas as pd
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from scipy.optimize import minimize

# Step 1. 加载数据
def loaddata(data):
     return np.loadtxt(data, delimiter=',')

data = loaddata('./data/sample_1.txt')
# 查看数据
print(data.shape)
print(data[:5, :])

# 运行结果：
(100, 3)
[[34.62365962 78.02469282  0.        ]
 [30.28671077 43.89499752  0.        ]
 [35.84740877 72.90219803  0.        ]
 [60.18259939 86.3085521   1.        ]
 [79.03273605 75.34437644  1.        ]]

# Step 2. 整理训练样本
# np.c_是按行连接两个矩阵，就是把两矩阵左右相加，要求行数相等
X = np.c_[np.ones((data.shape[0],1)), data[:,0:2]]
y = np.c_[data[:,2]]
print(X[:5])
print(X.shape)
print(y[:5])

# 运行结果：
[[ 1.         34.62365962 78.02469282]
 [ 1.         30.28671077 43.89499752]
 [ 1.         35.84740877 72.90219803]
 [ 1.         60.18259939 86.3085521 ]
 [ 1.         79.03273605 75.34437644]]
(100, 3)
[[0.]
 [0.]
 [0.]
 [1.]
 [1.]]

X的第一列是给偏置b预留的，其余两列为特征值。

# Step 3. 可视化
def plotData(data, label_x, label_y, label_pos, label_neg, axes=None):
    # 获得正负样本的下标(正样本表示考试通过的，负样本表示考试未通过)
    neg = data[:,2] == 0
    pos = data[:,2] == 1
    if axes == None:
        axes = plt.gca() # 获取图的坐标信息

    # 依次绘制正负样本，使用不同颜色区分
    axes.scatter(data[pos][:,0], data[pos][:,1], marker='*', c='r', s=30, linewidth=2, label=label_pos)
    axes.scatter(data[neg][:,0], data[neg][:,1], c='c', s=30, label=label_neg)
    axes.set_xlabel(label_x)
    axes.set_ylabel(label_y)
    axes.legend(frameon= True, fancybox = True);
    
plotData(data, 'English', 'Math', 'Pass', 'Fail')

从样本分布可以看出，正负样本分布的还算整齐，可以使用决策边界进行划分。

定义Sigmoid函数、损失函数（即目标函数）、梯度函数

# Step 4. 训练前准备
# Sigmoid函数
def sigmoid(z):
    return(1 / (1 + np.exp(-z)))

# 损失函数
def lossFunction(w, X, y):
    m = y.size
    h = sigmoid(X.dot(w.reshape(-1, 1)))
    # 套公式得到损失函数
    J = -1.0*(1.0/m)*(np.log(h).T.dot(y)+np.log(1-h).T.dot(1-y))
    if np.isnan(J[0]):
        return(np.inf)
    return J[0]

# 梯度计算函数
def gradient(X, y, w):
    grad = np.zeros(w.shape)
    error = (sigmoid(X.dot(w.reshape(-1, 1))).reshape(-1, 1) - y).ravel()
    for j in range(len(w.ravel())):
        # 循环更新w，得到梯度
        term = np.multiply(error, X[:,j])
        grad[j] = np.sum(term) / len(X)
    return grad

# 初始化w
initial_w = np.zeros(X.shape[1])
loss = lossFunction(initial_w, X, y)
grad = gradient(X, y, initial_w)
print('Loss: \n', loss)
print('Grad: \n', grad)

# 输出结果：
Loss: 
 [0.69314718]
Grad: 
 [ -0.1        -12.00921659 -11.26284221]

关于损失函数，代码里使用的是下面这个公式，而不是之前推导到最后的公式，因为用代码更好实现：

使用梯度下降法训练

# step 5. 使用梯度下降法训练
# 梯度下降函数
def gradientDescent(X, y, w, alpha=0.01, num_iters=1000):
    m = y.size
    J_history = np.zeros(num_iters)
    for iter in np.arange(num_iters):
        h = X.dot(w.reshape(-1, 1))
        grad = gradient(X, y, w)
        w = w - alpha*(1.0/m)*grad
        J_history[iter] = lossFunction(w, X, y)
    return(w, J_history)

# 画出每一次迭代和损失函数变化
w , Cost_J = gradientDescent(X, y, initial_w, alpha=0.15, num_iters=100000)

plt.plot(Cost_J)
plt.ylabel('Loss J')
plt.xlabel('Iter')

通过将学习率设置为0.15，并经过10W次的迭代后，损失函数计算出来的误差控制在了0.32左右。
下面对迭代出来的w进行可视化。

可视化决策边界

# Step 6. 可视化
plotData(data, 'English', 'Math', 'Pass', 'Fail')
# np.linspace(x1_min, x1_max) 生成的是x轴英语成绩的列向量
# np.linspace(x2_min, x2_max) 生成的是y轴数学成绩的列向量
# xx1, xx2为meshgrid生成的坐标矩阵
x1_min, x1_max = X[:,1].min(), X[:,1].max(),
x2_min, x2_max = X[:,2].min(), X[:,2].max(),
xx1, xx2 = np.meshgrid(np.linspace(x1_min, x1_max), np.linspace(x2_min, x2_max))
# 将预测值作为等高线绘制出来，即决策边界
# 将xx1和xx2作为x的特征值，得到预测值h
h = sigmoid(np.c_[np.ones((xx1.ravel().shape[0],1)), xx1.ravel(), xx2.ravel()].dot(w.T))
h = h.reshape(xx1.shape)
plt.contour(xx1, xx2, h, 1, linewidths=1, colors='b');

使用scipy.optimize.minimize实现

scipy.optimize.minimize可以根据现有的损失函数和梯度，优化整个梯度下降的过程，可以自动选择学习率，从而免去了不断手动调试的麻烦。

from scipy.optimize import minimize

res = minimize(lossFunction, initial_w, args=(X,y), jac=gradient, options={'maxiter':400})
print(res)

传入相应的参数，并设置迭代轮数（这里仅设置了400次），就可得到优化后的w。从输出可以看到，损失函数的误差已经减小到了0.2，比之前手动实现的要小很多，而且更加省时。

# Step 5. 可视化
plotData(data, 'English', 'Math', 'Pass', 'Fail')
x1_min, x1_max = X[:,1].min(), X[:,1].max(),
x2_min, x2_max = X[:,2].min(), X[:,2].max(),
xx1, xx2 = np.meshgrid(np.linspace(x1_min, x1_max), np.linspace(x2_min, x2_max))

h = sigmoid(np.c_[np.ones((xx1.ravel().shape[0],1)), xx1.ravel(), xx2.ravel()].dot(res.x.T))
h = h.reshape(xx1.shape)
plt.contour(xx1, xx2, h, 1, linewidths=1, colors='b');

从训练结果来看，手动和调包实现的效果好像差不多，但是明显后者更加省时省力。