回归GINI系数 回归系数sig

前言

之前对什么是逻辑回归，以及它的公式由来做了说明。在明确了该分类器的函数式后，那么最佳的回归系数是多少呢？这是值得思考的问题，本篇博客将会对这个问题进行探讨。

回顾逻辑回归公式

逻辑回归 = 线性回归 + sigmoid函数

线性回归：     z = w*x + b

sigmoid函数：y = ${1}\over{1+e^{-z}}$

逻辑回归：y = ${1}\over{1+e^{-w*x + b}}$

对于sigmoid函数，其输入z = $w_0$$x_0$ + $w_1$$x_1$ + $w_2$$x_2$ + … + $w_n$$x_n$，写成向量的形式为：z = $W^T$X，即将两个数值向量对应的元素相乘再求和得到z。这里的X就是样本数据，即分类器输入的数据；W就是需要求的最佳参数(系数)。分类器最终的分类效果好不好，精确率怎么样，都会受W的影响，所以要尽可能寻找最佳的W。寻找最佳W，会涉及到最优化理论的相关知识。

梯度上升法

这是一个最优化算法，这个算法基于的思想是：要找到某函数的最大值，最好的方法是沿着该函数的梯度方向探寻。如果梯度记作▽，那么函数f(x, y)的梯度有如下表示：
$\nabla f(x, y)=\left(\frac{\partial f(x, y)}{\partial x}, \frac{\partial f(x, y)}{\partial y}\right)$，这个梯度表示要沿着x的方向移动${\partial f(x, y)}\over{\partial x}$,沿着y的方向移动${\partial f(x, y)}\over{\partial y}$。需要注意的是，函数f(x, y)在待计算点必须有定义且可微。
梯度上升算法迭代公式：$w:=w+\alpha \nabla_{w} f(w)$
与之相反的是梯度下降算法，公式由加号变成减号：$w:=w-\alpha \nabla_{w} f(w)$，且用来求函数最小值。

python代码实现梯度上升找最佳参数

现在想要完成的任务是，在一些样本点上，利用梯度上升法找最佳回归系数，即找逻辑回归中最佳参数w。

1. 加载数据集

数据集：这是一个包含100个样本点的数据集，每个点包含两个特征（数值型）。即是一个100X3（100行3列）的矩阵，第三列是样本标签（类别）。

def load_dataset():
    """加载数据集"""
    dataSet = []  # 数据
    labelSet = []  # 标签
    fr = open('lr-testSet.txt')
    for line in fr.readlines():
        lineArr = line.strip().split()
        # 为了便于后续计算，将数据集第一列值设为1.0
        dataSet.append([1.0, float(lineArr[0]), float(lineArr[1])])
        labelSet.append(int(lineArr[2]))
    return dataSet, labelSet

2. sigmoid函数
   y = ${1}\over{1+e^{-z}}$

def sigmoid(X):
    """sigmoid函数"""
    return 1.0/(1+exp(-X))

3. 梯度上升法

这个函数有两个参数，第一个参数是100X3的矩阵，有100个样本，其中第一列${X_0}$是前面提到的便于后续计算，值设为1.0，后两列就是样本的2个特征${X_1}$和${X_2}$

def grad_ascent(dataSet, labelSet):
    """梯度上升算法"""
    dataMatrix = mat(dataSet)  # 转成numpy矩阵
    labelMatrix = mat(labelSet).transpose()  # 先转成numpy矩阵再进行矩阵转置
    m, n = shape(dataMatrix)
    print("样本数and特征数：", m, n)
    alpha = 0.001  # 向目标移动的步长
    maxCycles = 500  # 设置的最大的迭代次数
    weights = ones((n, 1))  # 表示回归系数，3*1的矩阵，值全为1
    for k in range(maxCycles):
        h = sigmoid(dataMatrix*weights)  # 矩阵相乘  # m*3 X 3*1
        error = (labelMatrix - h)  # m*1 - m*1
        # alpha * dataMatrix.transpose() * error表示在每一列上的一个误差情况
        # 得到x1,x2,xn的系数的偏移量
        weights = weights + alpha * dataMatrix.transpose() * error
    return weights

注意：ones函数可以用来创建任意维度和元素个数的数组，且元素值均为1；上面的代码中weights = ones((n, 1)) # 表示回归系数，3*1的矩阵，值全为1表示我们自行创建了一个n*1的矩阵(n是特征数，为3，其中后两个是特征，第一个是设置的$X_0$)。因为我们最终实际所求的是回归系数w，因为有两个特征再加上设置的$X_0$，所以我们最后得到的是一个3X1的矩阵，即是对应的回归系数。

4. 完整代码

from numpy import mat, shape, ones, exp  # 可以用from numpy import * 代替此句


def load_dataset():
    """加载数据集"""
    dataSet = []  # 数据
    labelSet = []  # 标签
    fr = open('lr-testSet.txt')
    for line in fr.readlines():
        lineArr = line.strip().split()
        if len(lineArr) == 1:
            continue
        # 为了便于后续计算，将数据集第一列值设为1.0
        dataSet.append([1.0, float(lineArr[0]), float(lineArr[1])])
        labelSet.append(int(lineArr[2]))
    return dataSet, labelSet


def sigmoid(inX):
    """sigmoid函数"""
    return 1.0 / (1 + exp(-inX))


def grad_ascent(dataSet, labelSet):
    """梯度上升算法"""
    dataMatrix = mat(dataSet)  # 转成numpy矩阵
    labelMatrix = mat(labelSet).transpose()  # 先转成numpy矩阵再进行矩阵转置
    m, n = shape(dataMatrix)
    print("样本数and特征数：", m, n)
    alpha = 0.001  # 向目标移动的步长
    maxCycles = 500  # 设置的最大的迭代次数
    weights = ones((n, 1))  # 表示回归系数，3*1的矩阵，值全为1
    for k in range(maxCycles):
        h = sigmoid(dataMatrix*weights)  # 矩阵相乘  # m*3 X 3*1
        error = (labelMatrix - h)  # m*1 - m*1
        # alpha * dataMatrix.transpose() * error表示在每一列上的一个误差情况
        # 得到x1,x2,xn的系数的偏移量
        weights = weights + alpha * dataMatrix.transpose() * error
    return weights


# 测试
dataSet, labelSet = load_dataset()
lr_w = grad_ascent(dataSet, labelSet)
print("最佳回归系数：\n", lr_w)

输出结果：

样本数and特征数： 100 3
最佳回归系数：
 [[ 4.12414349]
 [ 0.48007329]
 [-0.6168482 ]]

所以，现在求得的最佳的回归系数就是[ 0.48007329]和[-0.6168482 ]。

注意：
如果出现了如下错误信息：
return 1.0 / (1 + exp(-inX))
TypeError: only size-1 arrays can be converted to Python scalars
是因为在使用numpy数组的时候，我们最好也使用其中的数学函数处理数据。因为最开始在实现sigmoid函数的时候，exp我是从math模块中导入(from math import exp)的，就出现了上边的错误信息。修改方法就是从numpy模块中导入这个exp。