【机器学习】Logistic回归的梯度上升法

                       

一、问题引入
    对于目标值是连续变量的问题来说，使用线性回归可能会解决得很好，即使问题不能用线性模型描述时，也可以使用局部加权线性回归解决。但现实生活中有一种问题，输出值只有两种情况：yes or no.这类问题常见有：电子邮箱中的垃圾邮件分类(spam or not spam)，肿瘤为良性或者恶性等。在这些问题中，我们想预测的变量y，可以统一认为它只能取两个值0或1，这种问题叫分类(classification)问题，但只是最简单的二元分类问题，多类的问题待以后学习中讨论。如果我们尝试用线性回归来解决此类问题，碰巧的话，有时可能会解决得好，如图中粉红色线

但是如果出现了一个很远的样本点，训练出来的模型可能就变成蓝色线的样子。此时很明显分类的效果就很差。因此，应用线性回归来解决分类问题并不是一个好的想法。此时，我们便需要一种新的模型——logistic回归来解决分类问题。

二、问题分析
    对于输出值y∈{0, 1}的两类分类问题，我们作出一个假设：


函数g被称为logistic函数或sigmoid函数，至于为什么选择会选择这个函数，以后会涉及这个问题，暂时不深究。这个函数的图像是：

看起来有点像单位阶跃函数。根据这个函数，大于0.5的数据被划入1类，小于0.5的数据被归为0类。
有了这个函数，对于一个样本，我们可以得到它的概率分布：

综合起来就是：

此公式便为伯努利分布，这里的y∈{0, 1}.
现在我们就可以把问题转化为求logistic回归的最佳回归系数。由于logistic回归可以被看作是一种概率模型，且输出y发生的概率与回归参数θ有关，因此我们可以对θ进行最大似然估计(Maximum Likelihood Estimate)，使得y发生的概率最大，此时的θ便是最优的回归系数。对整个数据集求似然函数得：

为了计算方便，取似然函数的对数函数：

对上式运用梯度上升法，得到θ的迭代式：

求导过程不在此赘述，计算结果为：

此式便是本算法的核心公式！其中α为梯度上升步长，与梯度下降一样，决定了函数上升的快慢。

三、代码实现
1.Matlab版
(1) logistic回归的梯度上升算法，此处α设为0.001

function thetas = gradAscent(dataMat, classLabels, iterNum)    [m, n] = size(dataMat);    alpha = 0.001;    thetas = ones(n, 1)    for i = 1:1:iterNum        h = 1.0 ./ (1 + exp(-(dataMat * thetas)));      % sigmoid函数        error = classLabels - h;        thetas = thetas + alpha * dataMat' * error;     % 梯度上升    endend
  
  
1
2
3
4
5
6
7
8
9
10

(2)导入数据并绘图程序

clcclear all;data = load('testSet.txt');     % 导入数据dataMat = data(:,1:2); [m, n] = size(dataMat);dataSet = eye(m, n + 1);dataSet(:, 2:3) = dataMat;      % 输入x的样本dataSet(:, 1) = 1;              % 令第零维向量为1，便于计算labelMat = data(:,3);           % 输出y的样本j = 0;k = 0;xcord1 = [];ycord1 = [];xcord2 = [];ycord2 = [];for i = 1:1:m    if labelMat(i) == 1             % 对样本进行分类，归为1类        j = j + 1;        xcord1(j) = dataMat(i, 1);        ycord1(j) = dataMat(i, 2);    else                            % 归为0类        k = k + 1;        xcord2(k) = dataMat(i, 1);        ycord2(k) = dataMat(i, 2);    endendthetas = gradAscent(dataSet, labelMat, 500);x = -3.0 : 0.1 : 3.0;y = (-thetas(1) - thetas(2) * x) / thetas(3);     % 分类边界figure(1);xlabel('x1');ylabel('x2');plot(xcord1, ycord1, 'go', 'MarkerSize', 6, 'LineWidth', 2);hold on;plot(xcord2, ycord2, 'bx', 'MarkerSize', 6, 'LineWidth', 2);plot(x, y, '-r', 'LineWidth', 2)
  
  
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

(3)运行结果

可见，绿色点为1类，蓝色点为0类，红色直线是分类决策边界。分类效果比较理想，只有极个别点没有被分好。

2.Python版
(1) logistic回归的梯度上升算法，此处α设为0.001

def sigmoid(inx):    return 1.0 / (1 + exp(-inx));def gradAscent(dataMatIn, classLabels):    dataMatrix = mat(dataMatIn)    labelMat = mat(classLabels).transpose()    m, n = shape(dataMatrix)    alpha = 0.001    maxCycles = 500    theta = ones((n, 1))    for k in range(maxCycles):        h = sigmoid(dataMatrix * theta)        error = (labelMat - h)        theta = theta + alpha * dataMatrix.transpose() * error    return theta
  
  
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

(2)绘图程序如下：

def plotBestFit(theta):    dataMat,labelMat = ld.loadDataSet()    dataArr = array(dataMat)    n = shape(dataArr)[0]     xcord1 = []; ycord1 = []    xcord2 = []; ycord2 = []    for i in range(n):        if int(labelMat[i])== 1:            xcord1.append(dataArr[i,1]); ycord1.append(dataArr[i,2])        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='s')    ax.scatter(xcord2, ycord2, s = 30, c = 'blue')    x = arange(-3.0, 3.0, 0.1)    y = (-theta[0] - theta[1] * x) / theta[2]    ax.plot(x, y)    plt.xlabel('X1'); plt.ylabel('X2');    plt.show()
  
  
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

(3)运行结果

四、总结
    虽然算法实现起来只有短短十几行，但其中蕴含的数学知识是非常丰富的。从现在开始，之后的学习中会涉及到越来越多的概率论与数理统计的相关知识，幸好这学期有好好去学习这门课，最后的成绩也还不错。但是，课堂所传授的知识是远远不够，需要不断加以深入和扩展，才能使自己应对机器学习各种算法的推导和证明游刃有余。当然，线性代数也是很重要的工具，至少也要知道矩阵是怎么相乘的。由于我使用了两种语言去实现这个算法，在编码的过程中对两者的一些语法细节会有些混淆，此问题多写程序多加熟悉便可解决。两种语言实现起来并无太大差别，个人比较喜欢Python，可能是因为自己更偏向于当一个程序员，但是如果数学原理非常复杂，还是Matlab比较好使。总之，目前打算两种都要学好！