对数几率回归

对数几率回归

简要介绍

利用回归解决分类问题，也就是用一个函数拟合数据，最终将函数值映射到0或1 。搞懂对数几率回归，对理解神经网络很有帮助。

具体原理

给定数据集$D=(x_1,y_1),(x_2,y_2),\cdots (x_n,y_n).$其中每个样本输入有$n$个属性$(x_{11},x_{12},\cdots x_{1n})$，则会对应$n$个权重$(w_1,w_2,\cdots w_n)$，令$Z_1={\sum }_{i=1}^n{x}_{1i}{w}_i+b$，利用sigmod函数，最后分类结果为$\frac{1}{1+e^{-z}}$。

算法关键是如何选择参数，一个好的参数能使得对数损失函数
$$\begin{gathered} L(w,b)=-[y_{i}lo{g}_2\widehat{y_i}+(1-y_i)lo{g}_2(1-\widehat{y_i})]最小 \\ 其中y_i是第i个样本的类别，有0和1两种取值，\widehat{y_i}是预测分类. \end{gathered}$$
上面的损失函数其实只有$y_{i}lo{g}_2\widehat{y_i}$有效或者$(1-y_i)lo{g}_2(1-\widehat{y_i})$有效。当$y_i=\widehat{y_i}$时，损失函数为0最小，当$y_i与\widehat{y_i}$是0和1的取值时，即结果完全相反，此时损失函数为正无穷大。

损失函数为什么不直接选取最小二乘$\sigma (w,b)=|y_i-{\hat{y}}_i{|}^2$呢，因为最小二乘的损失函数不是凸的，而交叉熵$L(w,b)$是严格凸的，这样梯度下降可以到达最优解而不是局部解

如何使得损失函数最小，这里使用梯度下降的办法，梯度下降有随机梯度下降和平均梯度下降，下面考虑一个样本的梯度下降，有点类似随机梯度下降

给出各参数之间的关系：

$$\begin{gathered} L(w,b)=-[ylo{g}_2\hat{y}+(1-y)lo{g}_2(1-\hat{y})]. \\ \widehat{y_i}=\frac{1}{1+e^{-Z}}. \\ {Z}_i=\sum _{j=1}^n{x}_j{w}_j+b. \end{gathered}$$

从而求出梯度：

$$\begin{gathered} \frac{\partial L}{\partial w}=\frac{\partial L}{\partial \hat{y}}\ast \frac{\partial \hat{y}}{\partial Z}\ast \frac{\partial z}{\partial w} \\ =\frac{\hat{y}-y}{\hat{y}(1-y)}\ast \frac{1}{1+e^{-z}}(1-\frac{1}{1+e^{-z}})\ast x \\ =(\hat{y}-y)x. \\ 备注：此处w和x乃同维度的向量。 \\ 同理，\frac{\partial L}{\partial b}=\hat{y}-y. \end{gathered}$$

则可以更新参数：

$$\begin{gathered} w=w-\eta (\hat{y}-y)x. \\ b=b-\eta \hat{y}-y. \end{gathered}$$

代码实现

给出一部分数据格式

-0.017612 14.053064 0
-1.395634 4.662541 1
-0.752157 6.538620 0
-1.322371 7.152853 0
0.423363 11.054677 0
0.406704 7.067335 1
0.667394 12.741452 0
-2.460150 6.866805 1
0.569411 9.548755 0
-0.026632 10.427743 0

import numpy as npy
import matplotlib.pyplot as plt

# 读入数据函数
# 每一行前两个数字为x^1, x^2最后一个数据为标签y
# dataMat每一行为向量x = (1, x^1, x^2)
# labelMat只有一行，为对应的y
def loadDataSet():
    dataMat = []
    labelMat = []
    fr = open("testSet.txt")
    for line in fr.readlines():
        lineArr = line.strip().split()
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
        labelMat.append(int(lineArr[2]))
    return dataMat, labelMat

# 返回sigmoid函数计算结果
# inX为w*x
def sigmoid(inX):
    return 1.0 / (1 + npy.exp(-inX))

def gradAscent(dataMatIn, classLabels):
    dataMatrix = npy.mat(dataMatIn)
    # 将labelMat变成m*1的矩阵
    labelMat = npy.mat(classLabels).transpose()
    # m行n列
    m, n = npy.shape(dataMatrix)
    alpha = 0.001
    maxCycles = 500
    # 生成n*1列的矩阵
    weights = npy.ones((n, 1))
    for k in range(maxCycles):
        # m*n的矩阵乘上n*1的矩阵，结果为m*1的矩阵
        forecastLabel = sigmoid(dataMatrix * weights)
        # 两个m*1的矩阵相减
        error = (labelMat - forecastLabel)
        # n*m矩阵乘上m*1矩阵，得到n*1矩阵，在这里就是2*1
        # 这就是梯度上升法的alpha * L'(w)
        weights = weights + alpha * dataMatrix.transpose() * error
    return weights

# 该函数画出所有的数据点和最优的分割线
def plotBestFit(weights):
    dataMat, labelMat = loadDataSet()
    dataArr = npy.array(dataMat)
    n = npy.shape(dataArr)[0]
    xcord1 = []
    ycord1 = []
    xcord2 = []
    ycord2 = []
    for i in range(n):
        # 如果数据标签是y = 1
        if int(labelMat[i]) == 1:
            xcord1.append(dataArr[i, 1])
            ycord1.append(dataArr[i, 2])
        # 如果数据标签是y = 0
        else:
            xcord2.append(dataArr[i, 1])
            ycord2.append(dataArr[i, 2])
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(xcord1, ycord1, s = 30, c = "red", marker = "1")
    ax.scatter(xcord2, ycord2, s = 30, c = "green")

    # 画出分割线的图像即y = w * x
    # 原方程为0 = w0*x0 + w1*x1 + w2*x2
    x = npy.arange(-3.0, 3.0, 0.1)
    y = (-weights[0] - weights[1] * x) / weights[2]
    ax.plot(x, y)
    plt.xlabel("X1")
    plt.ylabel("X2")
    plt.show()

最终结果为：

[[ 4.12414349]
 [ 0.48007329]
 [-0.6168482 ]]

即$Z=4.12414349+0.48007329x_1-0.6168482x_2.$

决策结果为$sigmod(Z).$

---

以上内容参考《机器学习实战》Peter Harriongton著《机器学习》周志华著.