李宏毅机器学习logistic-regression

台大李宏毅机器学习课程第二次作业，logistic regression预测收入。这是我对这个作业的总结，方便日后复习。代码和训练数据已上传至github,戳我

知识准备

2019.11.24号这个周末用了两天的时间终于做完了，期间参考的各种资料我会在文章最后列出。那完成这次作业至少需要具备下面的知识点。

1. 微积分
2. 梯度下降原理
3. pandas和numpy的基本使用

作业要求

课程提供的作业要求文档.(其中的kaggle提交地址已失效,我在文章末尾提供了训练集和测试集的下载地址)。给定训练集和测试集，根据训练集建立模型判断测试集中的每个id对应的收入是否超过50K.

思路分析

首先，很明显的二分类问题，根据在b站(对，就是b站)课堂上学到内容选择使用logistic regression(对数几率回归)来实现.数值优化使用梯度下降，当然其他方法也是可以的。

1. 引入对数几率预测函数

先引入一个概念：广义线性模型。在西瓜书上的定义是这样的考虑单调可微函数$g(X)$令 $ y = g(w^Tx +b)$这样就得到了“广义线性模型(generalized linear model)“其中函数$g(X)$称为”联系函数(link function)“.显然，线性回归模型就是$g(X) = 1,x \in R$ 。 现在回到问题中，我们需要使用线性模型解决一个二分类的问题，那么只需要找到一个单调可微的函数$g(X)$将分类任务的真实标记y与线性回归模型的预测值联系起来即可。这句不理解没关系，接着看下一段就明白了

线性回归模型中产生的预测值$z = w^Tx + b$这里的$z$是一个实值$z \in R$。二分类问题中我们想要的预测值$y$是一个离散值0或者是1，即$y$的值是正例或者反例。那么我们只需要找一个函数 $ y = g(z) $ 它的输入值 $z$ 输出值是 $ y \in [0,1] $ 。最理想的函数是”单位阶跃函数“。如下：

$$g(z) = \begin{cases} 0 & z < 0\ 0.5&z=0 \ 1&z>0 \end{cases} $$

即若预测值$z$大于0判断为分类1，小于0则判断我分类2，预测值为临界值则可以任意判别。但是这个函数如下图它并不连续不能作为广义线性模型中的$g(X)$。于是我们希望找到一个函数一定程度上相似单位阶跃函数同时单调可微。对数几率函数就是这样一个函数：

$$g(z) = \frac{1}{1 + e^{-z}}$$

对数几率函数与单位阶跃函数图像如下：

我们把这个函数 代入到广义线性模型中就得到了对数几率回归的预测函数

$$y = \frac{1}{1+e^{-z}} = {1 \over 1+ e^{-{(w^Tx+b)}}}$$。

1.1 解释对数几率的含义

现在我们来解释一下“对数几率”这个词的含义。首先我们先把 $y = \frac{1}{1+e^{-z}} $ 中的$z$提取出来就得到了 $ln{\frac{y}{1-y}}=w^Tx+b$ 。其中如果把y看作为样本是正例的可能性，则$1-y$是其反例的可能性，这两者的比值 $ \frac{y}{1-y} $就称为几率，反映了$x$作为正例的可能性。对几率取对数得到的 $ ln\frac{y}{1-y} $就称为对数几率。

2.找到损失函数

现在已经知道了对数几率回归的预测函数，只需要找到损失函数，然后用梯度下降法确定w和b的值就完成任务了。对数几率回归中损失函数并是线性回归中的平方误差损失，而是交叉熵(可以参考这篇文章);或者也可以使用对数似然推导(可以参考书籍-机器学习-周志华第三章部分)出来损失函数结果是一样的。交叉熵老师在视频里已经给出了公式：

$$ L(w,b) = -[ylny^* + (1-y)ln(1-y^*)]$$

这个损失函数通常称作为 对数损失 (logloss)，这里的对数底为自然对数 ，其中真实值 是有 0/1 两种情况，而推测值$ y^* $ 由于借助对数几率函数，其输出是介于0~1之间连续概率值。仔细查看，不难发现，当真实值 时，第一项为0，当真实值 时，第二项为0，所以，这个损失函数其实在每次计算时永远都只有一项在发挥作用，那这不就可以转换为分段函数了吗，分段的形式如下:

不难发现，当真实值 为1时，输出值 越接近1，则 越小，当真实值 为 0 时，输出值 越接近于0，则 越小 (可自己手画一下 函数的曲线)。该分段函数整合之后就是上面我们所列出的 logloss 损失函数的形式。

3.使用梯度下降优化求解

我们确定了模型的损失函数，接下来就是根据这个损失函数，不断的优化模型的参数，从而获得拟合数据的最佳模型。损失函数如下：

$$L(w,b) = -[ylny^* + (1-y)ln(1-y^*)], 注意y不是自变量，y^*=\frac{1}{1+e^{-z}}, z = w^Tx + b$$

最小化损失值便得到

$(w^*,b^*) = argmin_{w,b}(L(w,b))$

梯度下降中的$(w,b)$的更新方式

$$w \leftarrow w - \alpha\frac{\partial L}{\partial w}$$

$$b \leftarrow b - \alpha\frac{\partial L}{\partial b}$$

由微积分的链式求导可以得到

$$\frac{\partial L}{\partial w} = \frac{\partial L}{\partial y^*} \frac{\partial y^*}{\partial z} \frac{\partial z}{\partial w} = (y^* - y)x$$

$$\frac{\partial L}{\partial b} = \frac{\partial L}{\partial y^*} \frac{\partial y^*}{\partial z} \frac{\partial z}{\partial b}=y^*-y$$

转化为矩阵运算：

为了简便运算我们在训练集上添加一个全为1的列，使偏置系数b更新的时候不受x值得影响，接着令$\beta = (W,b)$那么参数更新就可以写为

$$\beta = \beta - \alpha X^T(Y^* - Y)$$。

其中它们的shape如下：

$$\begin{cases} \beta & (n,1) \ X & (m,n) \ Y^*,Y&(m,1)\end{cases}$$

代码实现

下面是代码实现的部分，其中数据处理部分占了很大一部分，核心的梯度下降部分代码如下：

def gradient_secent(x_train, y_train):
    '''
    args:
        x_train: m x n
        y_train: m x 1
    return: wb

    '''
    lr_rate, threshold = 1, 0.01	# 初始化学习速度和阈值
    wb = np.ones((x_train.shape[1], 1))  # n x 1, 这里的wb就是刚刚的贝塔（就是那个和B很像的字母）
    g = np.zeros((x_train.shape[1], 1))	# 使用adagrad更新学习速率，
    y_predict = get_y_predict(wb, x_train)
    loss_value = get_loss_value(wb, y_train, y_predict)  # 计算损失值
    while True:
        wb, g = update_parameter(wb, lr_rate, x_train, y_train, y_predict, g)
        y_predict = get_y_predict(wb, x_train)
        pre_loss_value, loss_value = loss_value, get_loss_value(
            wb, y_train, y_predict)
        if abs(pre_loss_value - loss_value) < threshold: break
    return wb

完整的代码部分如下：

#!/usr/bin/env python
# encoding: utf-8
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

def binary_value(column_data,cond):
    '''
    args:
        column_data: 需要二值化的列
        cond: 二值化的条件，符合该条件的设置为1否则为0
    return: 返回处理好的数据
    '''
    column_data.mask(cond=cond,other=1,inplace=True)
    column_data.where(cond=column_data==1,other=0,inplace=True)
    #print(column_data)
    return column_data.astype('int64')

def one_hot_scale(data):
    '''
    对数据进行one-hot编码，并进行缩放.
    args:
        data: 原始数据
    return: 处理后的数据
    '''
    object_columns = [column for column in data.columns if data[column].dtype == 'object']
    number_columns = [column for column in data.columns if data[column].dtype == 'int64' and column != 'income' and column != 'sex']
    object_columns, number_columns = data[object_columns],data[number_columns]  # 提取相应的列
    
    # 对数据进行标准化缩放,缩放后值在0-1之间
    number_columns = (number_columns - number_columns.min()) / (number_columns.max() - number_columns.min())
    object_columns = pd.get_dummies(object_columns)
    data = pd.concat([object_columns, number_columns, data['sex'], data['income']],axis=1)
    
    return data

def data_pre_process(df_data,df_test_data):
    '''
    1. 进行数据预处理
    2. 返回训练集和验证集
    '''
    df_data,df_test_data = df_data.fillna(0),df_test_data.fillna(0)  # 对空值进行填充

    # 对性别和income列处理
    df_data['sex'] = binary_value(df_data['sex'].copy(deep=True),cond=df_data['sex']==' Male')
    df_test_data['sex'] = binary_value(df_test_data['sex'].copy(deep=True),cond=df_test_data['sex']==' Male')
    df_data['income'] = binary_value(df_data['income'].copy(deep=True),cond=df_data['income']==' >50K')

    df_data = one_hot_scale(df_data)
    df_test_data= one_hot_scale(df_test_data)

    # 把训练集分为训练集和验证集
    x_train = df_data.iloc[:int(len(df_data.index) * 0.80), :]
    x_valid = df_data.iloc[int(len(df_data.index) * 0.80):, :]
    
    y_train = np.array(x_train.iloc[:, -1:])
    x_train = np.array(x_train.iloc[:, :-1])
    x_train = np.c_[x_train, np.ones((x_train.shape[0], 1))]  # 添加一个横为1的列，用于简化偏置系数b的更新
    y_valid = np.array(x_valid.iloc[:, -1:])
    x_valid = np.array(x_valid.iloc[:, :-1])
    x_valid = np.c_[x_valid, np.ones((x_valid.shape[0], 1))]
    
    y_test = np.array(df_test_data.iloc[:,-1:])
    x_test = np.array(df_test_data.iloc[:,:-1])
    x_test = np.c_[x_test, np.ones((x_test.shape[0], 1))]  # 添加一个横为1的列，用于简化偏置系数b的更新
    return x_train, y_train, x_valid, y_valid,x_test,y_test


def gradient_secent(x_train, y_train):
    '''
    args:
        x_train: m x n
        y_train: m x 1
    return: wb
    '''
    lr_rate, threshold = 1, 5
    wb = np.ones((x_train.shape[1], 1))  # n x 1
    g = np.zeros((x_train.shape[1], 1))
    y_predict = get_y_predict(wb, x_train)
    loss_value = get_loss_value(wb, y_train, y_predict)  # 计算损失值
    loss_value_list = []
    while True:
        loss_value_list.append(loss_value)
        wb, g = update_parameter(wb, lr_rate, x_train, y_train, y_predict, g)
        y_predict = get_y_predict(wb, x_train)
        pre_loss_value, loss_value = loss_value, get_loss_value(wb, y_train, y_predict)
        if abs(pre_loss_value - loss_value) < threshold: break  # 循环结束条件
    plot_loss_value(loss_value_list)    #画图
    return wb


def update_parameter(wb, lr_rate, x_train, y_train, y_predict, g):
    '''
    args:
        wb: n x 1
        g: adagrad 里的 G, n x 1
        lr_rate: 学习速率
        y_train: 真实值，m x 1
        y_predict: 预测值，m x 1
    return: wb, g
    '''
    descent = x_train.T.dot(y_predict - y_train)  # n x 1
    g = g + descent * descent
    lr_rates = lr_rate / (np.sqrt(g + 0.00000001))  # n x 1
    wb = wb - lr_rates * descent
    return wb, g


def get_y_predict(wb, x_train):
    '''
    args:
        wb: 参数，n x 1
        x_train: 训练集，m x n
    return:计算预测值 y.shape = m x 1
    '''
    y = wb.T.dot(x_train.T)  # 1 x n * n x m = 1 x m
    y = 1 / (np.exp(-1 * y) + 1)  # 1 x m
    y = y.T  # m x 1
    return y

def get_loss_value(wb, y_train, y_predict):
    '''
    args:
        wb: 参数，n x 1
        y_train: 真实值，m x 1
        y_predict: 预测值，m x 1
    return: loss_value, 损失值
    '''
    res = y_train.T.dot(np.log(y_predict))  # 1 x m * m x 1
    res = -1 * res
    tmp = (y_train - 1).T.dot(np.log(1 - y_predict))
    res = res + tmp
    return res[0][0]


def get_acc(wb, x, y):
    '''
    args:
        wb: n x 1
        x: m x n
        y: m x 1
    return: 正确率 acc
    '''
    y_predict = get_y_predict(wb, x)  # m x 1
    y_predict[y_predict > 0.5] = 1
    y_predict[y_predict <= 0.5] = 0
    acc = sum(y_predict == y)[0] / y.shape[0]
    return acc

def test_data_process(df_data):
    '''
    args:
        df_data: m x n
    return:
        x_test, y_test
    '''

def main():
    train_data = pd.read_csv('./data/train.csv')
    test_data = pd.read_csv('./data/test.csv')
    x_train, y_train, x_valid,y_valid,x_test,y_test = data_pre_process(train_data,test_data)
    wb = gradient_secent(x_train, y_train)
    valid_acc = get_acc(wb, x_valid, y_valid)
    print("valid accuracy:{}".format(valid_acc))
    test_acc = get_acc(wb,x_test,y_test)
    print("tet accuracy:{}".format(test_acc))
    plt.show()

def plot_loss_value(loss_values):
    '''
    画出loss value的变化图
    '''
    plt.plot(range(len(loss_values)),loss_values,linewidth=4,color='g')

if __name__ == '__main__':
    '''
    李宏毅机器学习第二次作业
    1. 对数几率回归
    2. 使用梯度下降寻找最优解
    '''
    print("正在执行中，可能需要耗费大量时间，如果超过10分钟未响应请手动关掉即可。")
    main()

运行截图：

可以看出在验证集和测试集中正确率还算看的过去。。。可能做了特征工程后正确率会高一点。

参考链接

感谢下面的博主和作者。

1. 《机器学习》——周志华

2. 对数几率回归 —— Logistic Regression

3. 数据集连接