《机器学习实战》第5章逻辑回归算法的C++实现方式

在VS2013平台下，调用armadillo线性代数库实现《机器学习实战》第5章的算法。算法的原理和数学推导请参考书籍或相关博客http://blog.youkuaiyun.com/zouxy09/article/details/20319673。
由于在C++控制台上不易画图，这里将训练好的权值保存到.txt文档中，在保存目录下建立一个.m脚本用matlab展示可视化效果。

C++代码

这里函数的功能和命名均与《机器学习实战》中的函数保持高度一致

#include "stdafx.h"
#include "armadillo"
#include <iostream>
using namespace arma;
using namespace std;

//加载数据
bool loadDataSet(mat &dataMat,vec & labelVec)
{
    bool openSuccess = true;
    dataMat = ones<mat>(500, 3);
    labelVec = zeros<vec>(500);
    fstream in("D:\\TestCode\\Machine-Learning-In-Action\\Chapter5\\testSet.txt");
    if (!in.is_open())
    {
        openSuccess = false;
        return openSuccess;
    }
    int dataNum = 0; //数字总数
    while (!in.eof())
    {
        dataNum++;
        int rowNum = int(dataNum / 3);//应分配的行,从0行开始
        int colNum = int(dataNum % 3);//应分配的列
        if (colNum != 0)
        {
            in >> dataMat(rowNum, colNum);
        }
        else
        {
            in >> labelVec(rowNum - 1);//(分析：当dataNum=3时rowNum=1，但此时应为第0行的标签)
        }
    }
   //截取矩阵中有数据部分
    int linesNum = int(dataNum / 3);
    dataMat = dataMat.rows(0, linesNum - 1);
    labelVec = labelVec.rows(0, linesNum - 1);
    return openSuccess;
}

//sigmoid行数，输入和返回值都是列向量
vec sigmoid(vec inX)
{
    return 1.0 / (1.0 +exp(-inX));
}
//梯度上升函数（批处理）
vec gradAscent(mat dataMatIn, vec classLabels)
{
    double alpha = 0.001;//梯度上升步长
    int maxCycles = 500;//最大循环次数
    int cols = dataMatIn.n_cols;//每个样本的特征数
    vec weights = ones<vec>(cols);//形成初始权值，全为1
    for (size_t i = 0; i < maxCycles; i++)
    {
        vec h = sigmoid(dataMatIn*weights);
        vec error = (classLabels - h);
        weights = weights + alpha * dataMatIn.t()*error;
    }
    return weights;
}
//随机梯度上升算法
vec stocGradAscent(mat dataMatIn, vec classLabels)
{
    double alpha = 0.01;
    int cols = dataMatIn.n_cols;//每个样本的特征数
    int rows = dataMatIn.n_rows;//样本个数
    vec weights = ones(cols);//初始权值
    for (size_t i = 0; i < rows; i++)
    {
        vec h = sigmoid(dataMatIn.row(i)*weights);//起始只有一个值，标量
        vec error = classLabels(i) - h;
        weights = weights + alpha*dataMatIn.row(i).t()*error;
    }
    return weights;
}

//改进的随机梯度上升算法
vec stocGradAscent1(mat dataMatIn, vec classLables, int numIter = 150)
{
    int cols = dataMatIn.n_cols;//每个样本的特征数
    int rows = dataMatIn.n_rows;//样本个数
    vec weights = ones(cols);//初始权值
    for (int j = 0; j < numIter; j++)
    {
        uvec dataIndex = linspace<uvec>(0, rows - 1, rows);//产生0到rows-1的序列
        for (int i = 0; i < rows; i++)
        {
            double alpha = 4.0 / (1.0 + i + j) + 0.01;//改进处1：逐渐减小的步长
            int randIndex = rand() % (dataIndex.n_elem);//产生一个随机数
            vec h = sigmoid(dataMatIn.row(randIndex)*weights);//起始只有一个值，标量
            double error = classLables(randIndex) - h(0);
            weights = weights + alpha*dataMatIn.row(randIndex).t()*error;
            dataIndex.shed_row(randIndex);
        }
    }
    return weights;
}


int _tmain(int argc, _TCHAR* argv[])
{
    mat dataSet;
    vec labelVec;
    bool success = loadDataSet(dataSet, labelVec);
    //vec weight = gradAscent(dataSet, labelVec);//批处理梯度上升
    //vec weight = stocGradAscent(dataSet, labelVec);//随机梯度上升算法
    vec weight = stocGradAscent1(dataSet, labelVec);//改进的随机梯度上升算法
    weight.print("weight");
    //将训练好的权值存到weight.txt中，以供matlab绘图
    fstream out("D:\\TestCode\\Machine-Learning-In-Action\\Chapter5\\weight.txt",ios_base::out);
    for (int i = 0; i < weight.n_rows; i++)
    {
        out << weight(i)<<" ";
    }
    uvec a = linspace<uvec>(0, 4, 5);

    system("pause");
    return 0;
}

---

Matlab绘图脚本

%% 从txt中读取数据并绘图
testSet = importdata('testSet.txt');
weight = importdata('weight.txt');
testSet_x0 = testSet((testSet(:,3) < 1),1); %提取标签为0的x
testSet_y0 = testSet((testSet(:,3) < 1),2); %提取标签为0的y

testSet_x1 = testSet((testSet(:,3) > 0),1); %提取标签为1的x
testSet_y1 = testSet((testSet(:,3) > 0),2); %提取标签为1的y
%% 绘图
plot(testSet_x0,testSet_y0,'bo');
hold on;
plot(testSet_x1,testSet_y1,'r*');
hold on;
Xrange= minmax(testSet(:,1)');% x的分布范围
x = Xrange(1)-0.5 : 0.2 : Xrange(2)+0.5;
y = (-weight(1,1)-weight(1,2)*x)/weight(1,3);
plot(x,y);
xlabel('x');
ylabel('y');
legend('0','1');

算法效果

1.调用批处理梯度上升

2.调用随机梯度上升

3.调用改进的随机梯度上升