Logistic回归与牛顿法（附Matlab实现）

回归，是一种连续模型，受噪声的影响较大，一般都是用来做预测的，但也有除外，比如本文要讲的Logistic回归就是用来做分类的。

Logistic Regress

Logistic一般用于二分类问题，不同于之前讲的线性回归，它是用一条直线来分割两种不同类别的样本。其函数形式为



若令

，则可以等价为



这个函数也叫sigmoid函数。其函数图像如下



得到了Logistic回归模型后，接着就要找到合适的θ去拟合它了。

首先假设



可得概率公式



再假设训练集独立同分布，类似于线性回归，这里也用最大似然估计得



将其对数化



接着求导



其中上式第一行到第二行之间其实有一个很有意思的东西，g（θ）的一个性质



之后也和雷同的和上文一样，这里可以直接用梯度下降法来求解



几乎一模一样，只是其中的h（x）变了，其实Logistic回归和线性回归都是属于同一个模型-----广义线性模型（Generalized Linear Models，GLM），这个具体下一次再写。

Newton's Method

除了用梯度下降，我们还可有用一个新的方法来解决这个问题，一个更快的方法----牛顿法（Newton’s method）

先介绍下这种方法，牛顿法主要是找到一个函数等于0的值，即当f（θ）=0时，θ的取值。它的更新公式为



至于它的具体操作用文字很难说清楚（文笔太烂了），这接用例子说明吧



从最左图看起，其中直线y=0和函数f（θ）相交于θ=1.3的位置，其中1.3就是我们要求的东西。这里假设θ从初始值为4.5开始更新（如中间图），作θ=4.5在函数f（θ）上的切线，这条切线和y=0的交点处的θ值就是更新后的θ位置（从最右图可以看出），暂且将更新后的θ定义为θ1，从图中可得f’（θ）=f（θ）/（θ - θ1），进而可得（θ - θ1）= f（θ）/f’（θ），很眼熟，和公式中的最后一项一样，这里f（θ）/f’（θ）就是每一次要更新的差值，反复迭代，直到新的θ等于0。

说完方法，就得开始用该方法解决问题了，这里用牛顿法来求最大似然函数l（θ）的极大值，即l’（θ）=0。公式变为



又由于θ是一个向量值，最后的更新表达式为



其中H为hessian矩阵，其形式为



PS：与梯度下降法相比，牛顿法的迭代次数更少，收敛的更快，但其每一次的花费却比梯度下降发更大，因为其中要处理hessian矩阵。

In Matlab

下面是用Matlab实现牛顿法的具体代码，数据集是乱找的，在代码中就不显示出来了。

% Logstic Regress with Newton's Method

clear all; close all; clc

% data
% x(20*2)为数据的特征矩阵
x = [0.230000 0.394000;
.......
0.780000 0.282000];
% y（20*1）为数据的标记矩阵
y = [0;
......
1];

[m, n] = size(x);

% x0 = 1
x = [ones(m, 1), x];

% plot the datas
figure
pos = find(y); neg = find(y == 0);%find是找到的一个向量，其结果是find函数括号值为真时的值的编号
plot(x(pos, 2), x(pos,3), '+')
hold on
plot(x(neg, 2), x(neg, 3), 'o')
hold on
xlabel('axis X')
ylabel('axis Y')

% Initialize fitting parameters
theta_new = zeros(n+1, 1);
theta_old = ones(n+1, 1);
e = 0.00000000001;

% Define the sigmoid function
g = inline('1.0 ./ (1.0 + exp(-z))');

% Newton's method
while ((theta_new - theta_old)*(theta_new - theta_old)' > e)
theta_old = theta_new;
% Calculate the hypothesis function
z = x * theta_new;
h = g(z);

% Calculate gradient and hessian.
grad = x' * (h-y);
H = x' * diag(h) * diag(1-h) * x;

theta_new = theta_new - H\grad;
end

theta_new

plot_x = [x(:,2)];
plot_y = (-1./theta_new(3)).*(theta_new(2).*plot_x +theta_new(1));
plot(plot_x, plot_y)
legend('Label 1', 'Label 2', 'Decision Boundary')
hold off

下图是未经处理的数据集plot出来的图。



下图是处理后的图



思考

尽管Logistic回归和线性回归的目的并不相同，一个是分类，一个是预测，但从本质上讲Logistic回归仍然是一种线性回归，只是在回归的结果上又加上了一层映射，将其映射到可分类的情况（其中带入到sigmoid函数的是

，而z在线性回归中是等于h（x）的）。g(z)可以将连续值映射到0和1上。

疑问

Hessian矩阵的向量表达式是怎么来的？