多维牛顿法matlab,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矩阵的向量表达式是怎么来的？