机器学习基础2-二项逻辑回归及C++实现

分类问题

生活中我们会遇到很多分类问题,如一封邮件是否是垃圾邮件等。那么能否建立一种模型来给出垃圾邮件的概率?通过前面线性回归一文我们知道线性回归的结果是一个连续值,并且它的范围是无法限定的,而分类问题需要我们给出一个0.0~1.0之间的概率值,那么有没有什么办法可以将线性回归的结果值映射为一个0.0~1.0之间的概率值呢?有!它就是逻辑函数。

逻辑函数是一种S函数,函数形式为:

$$ p(t)=\frac{1}{1+e^{-t}} $$

逻辑函数的函数图像为(来源于维基百科):

二项逻辑回归模型

利用逻辑函数我们可以构造二项逻辑回归模型(binomial logistic regression model)。二项逻辑回归模型中的二项意为二分类,也就是分类类别 $y\in{0,1}$ 。

设特征数量为 $n$ 即,特征为 $x$ ,特征的参数为 $\theta$ ,我们定义如下:

\begin{align} \\ x & = \begin{bmatrix} 1 \ x_1 \ \cdots \ x_n \end{bmatrix} \\ \\ \theta & = \begin{bmatrix} \theta_0 \\ \theta_1 \\ \vdots \\ \theta_n \end{bmatrix} \\ \\ x\theta & = \theta_0 + \theta_1x_1 + \theta_2x_2 + \cdots + \theta_nx_n \\ \\ h_\theta(x) & = \frac{1}{1+e^{x\theta}} \end{align}

$h_\theta(x)$ 就是我们的概率函数,我们可以假设它为分类类别为1的概率,那么类别1和类别0的概率分别为:

\begin{align} \\ P(y=1 \mid x) & =h_\theta(x) \\ \\ P(y=0 \mid x) & =1-h_\theta(x) \\ \end{align}

综合起来可以写成: \begin{align} \\ P(y \mid x) & =h_\theta(x)^y(1-h_\theta(x))^{1-y} \\ \end{align}

设训练样本数为 $m$,训练样本集为 $X$ ,训练输出集为 $Y$ ,如下: \begin{align} X & = \begin{bmatrix} x^{0} \\ x^{1} \\ \cdots \\ x^{m-1} \end{bmatrix} \\ \\ Y & = \begin{bmatrix} y^{0} \\ y^{1} \\ \vdots \\ y^{m-1} \end{bmatrix} \\ \end{align}

我们的目标是已知 $X$ 和 $Y$ 的情况下得到最优的 $\theta$。

似然函数

哪个 $\theta$ 是最优的?我们需要先定义似然函数:

\begin{align} \\ L(\theta) &= \prod_{i=0}^{m-1} P(y^i \mid x^i) \\ \\ L(\theta) &= \prod_{i=0}^{m-1}h_\theta(x^i)^{y^i} (1-h_\theta(x^i))^{1-y^i} \\ \end{align}

我们在似然函数中引入自然对数以方便后续的求导,则:

\begin{align} \\ L(\theta) &= \log(\prod_{i=0}^{m-1}h_\theta(x^i)^{y^i} (1-h_\theta(x^i))^{1-y^i}) \\ L(\theta) &= \sum_{i=0}^{m-1}\log(h_\theta(x^i)^{y^i} (1-h_\theta(x^i))^{1-y^i}) \\ L(\theta) &= \sum_{i=0}^{m-1}(y^i\log(h_\theta(x^i)) + (1-y^i)\log(1-h_\theta(x^i))) \\ \end{align}

很明显似然函数最大值对应的 $\theta$ 就是我们求解的目标,所以问题变为: $$ \max_\theta L_\theta $$

梯度上升法

使用梯度上升法可以帮助我们找到似然函数的最大值,参数 $\theta_j$的梯度为:

$$ \frac{\partial L(\theta)}{\partial \theta_j}=\frac{\partial}{\partial \theta_j} \sum_{i=0}^{m-1}(y^i\log(h_\theta(x^i)) + (1-y^i)\log(1-h_\theta(x^i))) $$

假设样本数为1,则得到如下: \begin{align} \\ \frac{\partial L(\theta)}{\partial \theta_j} &= \frac{\partial}{\partial \theta_j} (y\log(h_\theta(x)) + (1-y)\log(1-h_\theta(x))) \\ \\ \frac{\partial L(\theta)}{\partial \theta_j} &= \frac{y}{h_\theta(x)} \frac{\partial}{\partial \theta_j} h_\theta(x) + \frac{1-y}{1-h_\theta(x)}\frac{\partial}{\partial \theta_j}(1-h_\theta(x)) \\ \\ \frac{\partial L(\theta)}{\partial \theta_j} &= (\frac{y}{h_\theta(x)} - \frac{1-y}{1-h_\theta(x)})\frac{\partial}{\partial \theta_j} h_\theta(x) \\ \end{align}

又因为:

\begin{align} \\ \frac{\partial}{\partial \theta_j} h_\theta(x) &= \frac{\partial}{\partial \theta_j} \frac{1}{1+e^{x\theta}} \\ \\ \frac{\partial}{\partial \theta_j} h_\theta(x) &=\frac{-1}{(1+e^{x\theta})^2} \frac{\partial}{\partial \theta_j}e^{x\theta} \\ \\ \frac{\partial}{\partial \theta_j} h_\theta(x) &=\frac{-e^{x\theta}}{(1+e^{x\theta})^2} \frac{\partial}{\partial \theta_j} x\theta \\ \\ \frac{\partial}{\partial \theta_j} h_\theta(x) &=-h_\theta(x)(1-h_\theta(x)) \frac{\partial}{\partial \theta_j} x\theta \\ \\ \frac{\partial}{\partial \theta_j} h_\theta(x) &=-h_\theta(x)(1-h_\theta(x)) x_j \\ \end{align}

所以:

\begin{align} \\ \frac{\partial L(\theta)}{\partial \theta_j} &= (\frac{y}{h_\theta(x)} - \frac{1-y}{1-h_\theta(x)})(-h_\theta(x)(1-h_\theta(x)) x_j) \\ \\ \frac{\partial L(\theta)}{\partial \theta_j} &= -(y(1-h_\theta(x)) - (1-y)h_\theta(x))x_j \\ \\ \frac{\partial L(\theta)}{\partial \theta_j} &= -(y-h_\theta(x))x_j \\ \end{align}

多个样本的正确公式为: \begin{align} \\ \frac{\partial L(\theta)}{\partial \theta_j} &= -\sum_{i=0}^{m-1} (y^i-h_\theta(x^i))x_j^i \\ \end{align}

梯度上升法下的 $\theta_j$ 的更新公式为:($\alpha$ 为学习速度)

\begin{align} \\ \theta_j &:= \theta_j + \alpha \frac{\partial L(\theta)}{\partial \theta_j} \\ \\ \theta_j &:= \theta_j - \alpha \sum_{i=0}^{m-1} (y^i-h_\theta(x^i))x_j^i \\ \end{align}

关于随机梯度,批量梯度以及学习速度在上一文线性回归中我们已经介绍过,所以这边不再解释。

C++代码实现

我们定义如下的接口:

    typedef LMatrix<float> LRegressionMatrix;

    /// @brief 回归中的ZERO分类
    #ifndef REGRESSION_ZERO
    #define REGRESSION_ZERO 0.0f
    #endif

    /// @brief 回归中的ONE分类
    #ifndef REGRESSION_ONE
    #define REGRESSION_ONE 1.0f
    #endif

    /// @brief 逻辑回归(分类)
    class LLogisticRegression
    {
    public:
        /// @brief 构造函数
        LLogisticRegression();

        /// @brief 析构函数
        ~LLogisticRegression();

        /// @brief 训练模型
        /// 如果一次训练的样本数量为1, 则为随机梯度下降
        /// 如果一次训练的样本数量为M(样本总数), 则为梯度下降
        /// 如果一次训练的样本数量为m(1 < m < M), 则为批量梯度下降
        /// @param[in] xMatrix 样本矩阵, 每一行代表一个样本, 每一列代表样本的一个特征
        /// @param[in] yVector(列向量) 样本标记向量, 每一行代表一个样本, 值只能为REGRESSION_ONE或REGRESSION_ZERO 
        /// @param[in] alpha 学习速度, 该值必须大于0.0f
        /// @return 成功返回true, 失败返回false(参数错误的情况下会返回失败)
        bool TrainModel(IN const LRegressionMatrix& xMatrix, IN const LRegressionMatrix& yVector, IN float alpha);

        /// @brief 使用训练好的模型预测数据
        /// @param[in] xMatrix 需要预测的样本矩阵
        /// @param[out] yVector 存储预测的结果向量(列向量), 值为REGRESSION_ONE标记的概率
        /// @return 成功返回true, 失败返回false(模型未训练或参数错误的情况下会返回失败)
        bool Predict(IN const LRegressionMatrix& xMatrix, OUT LRegressionMatrix& yVector) const;

        /// @brief 计算似然值, 似然值为0.0~1.0之间的数, 似然值值越大模型越好
        /// @param[in] xMatrix 样本矩阵, 每一行代表一个样本, 每一列代表样本的一个特征
        /// @param[in] yVector(列向量) 样本输出向量, 每一行代表一个样本, 值只能为REGRESSION_ONE或REGRESSION_ZERO
        /// @return 成功返回似然值, 失败返回-1.0f(参数错误的情况下会返回失败)
        float LikelihoodValue(IN const LRegressionMatrix& xMatrix, IN const LRegressionMatrix& yVector) const;

    private:
        CLogisticRegression* m_pLogisticRegression; ///< 逻辑回归实现类
    };

LMatrix是我们自定义的矩阵类,用于方便机器学习的一些矩阵计算,关于它的详细代码可以查看链接:猛戳我

我们为LLogisticRegression设计了三个方法TrainModel,Predict以及LikelihoodValue,用于训练模型,预测新数据以及计算似然值。

我们看一下TrainModel的实现:

    LRegressionMatrix X;
    Regression::SamplexAddConstant(xMatrix, X);

    const LRegressionMatrix& Y = yVector;

    LRegressionMatrix& W = m_wVector;
    LRegressionMatrix XT = X.T();

    /*
    如果h(x)  =  1/(1 + e^(X * W)) 则
    wj = wj - α * ∑((y - h(x)) * xj)
    如果h(x)  =  1/(1 + e^(-X * W)) 则
    wj = wj + α * ∑((y - h(x)) * xj)
    */

    LRegressionMatrix XW(X.RowLen, 1);
    LRegressionMatrix DW;

    LRegressionMatrix::MUL(X, W, XW);
    for (unsigned int m = 0; m < XW.RowLen; m++)
    {
        this->Sigmoid(XW[m][0], XW[m][0]);
    }

    LRegressionMatrix::SUB(Y, XW, XW);
    LRegressionMatrix::MUL(XT, XW, DW);

    LRegressionMatrix::SCALARMUL(DW, -1.0f * alpha, DW);
    LRegressionMatrix::ADD(W, DW, W);

逻辑函数定义如下:

    /// @brief S型函数
    /// @param[in] input 输入值
    /// @param[out] output 存储输出值
    void Sigmoid(float input, float& output) const
    {
        output = 1.0f/(1.0f + exp(input));
    }

我们看一下Predict的实现:

    LRegressionMatrix X;
    Regression::SamplexAddConstant(xMatrix, X);

    yVector.Reset(X.RowLen, 1, 0.0f);
    LRegressionMatrix::MUL(X, m_wVector, yVector);

    for (unsigned int m = 0; m < yVector.RowLen; m++)
    {
        this->Sigmoid(yVector[m][0], yVector[m][0]);
    }

最后就是LikelihoodValue的实现:

    LRegressionMatrix predictY;
    this->Predict(xMatrix, predictY);

    float likelihood = 1.0f;
    for (unsigned int i = 0; i < yVector.RowLen; i++)
    {
        if (yVector[i][0] == REGRESSION_ONE)
            likelihood *= predictY[i][0];
        else if (yVector[i][0] == REGRESSION_ZERO)
            likelihood *= (1.0f - predictY[i][0]);
        else
            return -1.0f;
    }

    return likelihood;

以上完整的代码可以在链接:猛戳我查看,我们的逻辑回归被定义在文件LRegression.h和LRegression.cpp中。

其他章节链接

机器学习基础1-线性回归及C++实现

机器学习基础2-二项逻辑回归及C++实现

机器学习基础3-Softmax回归及C++实现


赞助作者写出更好文章


您还未登录,登录GitHub账号发表精彩评论

 GitHub登录


最新评论

  • strutsjava4560

    发表于: 2017年12月19日  

    回复@burnell liu

    公式排版有误?

  • burnell liu (作者)

    发表于: 2017年12月21日  

    回复@strutsjava4560

    求教哪里有误

 

 

刘杰

28岁, 现居苏州

微信:

BurnellLIU

邮箱:

burnell_liu@outlook.com

Github:

https://github.com/BurnellLiu

简介:

努力做一个快乐的程序员, good good study, day day up!

本站: 登录 注册   友情链接: Will Mao

苏ICP备16059872号. Copyright © 2017. http://www.coderjie.com. All rights reserved.

账号登录

注册 忘记密码