📚描述

在本练习中，我们将使用牛顿法对分类问题实现逻辑回归。

📚数据

在本练习中，假设一所高中有一个数据集，该数据集表示40名被大学录取的学生和40名未被录取的学生。每个 $x^{ ( i )} , y^{ ( i )} )$ 训练示例包含学生在两个标准化考试中的成绩以及学生是否被录取的标签。我们的任务是建立一个二元分类模型，根据学生在两门考试中的成绩来估计大学录取的机会。在训练数据中

x 数组的第一列代表所有测试1的分数，第二列代表所有测试2的分数。
y 矢量使用" 1“标记被录取的学生和" 0”标记没被录取的学生。

📚绘制数据

将训练样例的数据加载到程序中，并将 $x_0 = 1$ 截距项添加到x 矩阵中。在开始牛顿方法之前，我们将首先使用不同的符号绘制数据来表示两个类。在Matlab / Octave中，我们可以使用find命令分离正类和负类：

% find返回满足指定条件的行的索引
pos = find ( y == 1 ) ; neg = find ( y == 0 ) ;
% Assume the features are in the 2nd and 3rd columns of x
plot ( x ( pos , 2 ) , x ( pos , 3 ) , '+' ) ; hold on
plot ( x ( neg , 2 ) , x ( neg , 3 ) , ' o ' )

x = load('ex4x.dat');
y = load('ex4y.dat');
m = length(y);
x = [ones(m,1),x] ;

% find返回满足指定条件的行的索引
pos = find ( y == 1 ) ; neg = find ( y == 0 ) ;
% Assume the features are in the 2nd and 3rd columns of x
plot ( x ( pos , 2 ) , x ( pos , 3 ) , '+' ) ; hold on
plot ( x ( neg , 2 ) , x ( neg , 3 ) , ' o ' )
xlabel('Exam1 score');
ylabel('Exam2 score');
legend('Admitted','Not admitted');

在这里插入图片描述

📚牛顿法

回想一下逻辑回归，假设函数为

在这里插入图片描述

Matlab/Octave没有sigmoid函数库，所以你必须自己定义它。最简单的方法是通过内联表达式。

g = inline('1.0 ./ (1.0 + exp(-z))');
%Usage: To find the value of sigmoid
%evaluated at 2, call g(2)

代价函数为

在这里插入图片描述

我们的目标是使用牛顿的方法来最小化这个函数。回想一下，牛顿方法的更新规则是

在这里插入图片描述

在逻辑回归中，梯度和海森值是

在这里插入图片描述

注意，上面给出的公式是向量化的版本。具体来说，这意味着 $x^{(i)}∈R^{n+1}$ , $x^{(i)}(x^{(i)})^T∈R^{(n+1)×(n+1)}$ 其中 $h_θ ( x ^{( i )} )$ 和 $y^{( i )}$ 是标量。

现在，在你的程序中实现牛顿法，从初值θ = 0 开始。为了确定使用多少迭代，计算每次迭代过程中的 $L (θ)$ 和绘制你的结果。牛顿方法通常在5-15次迭代中收敛。如果您发现自己使用了更多的迭代，那么您应该检查实现中的错误。

收敛后，使用theta值来找到分类问题中的决策边界，它对应于 $θ^Tx=0$

在这里插入图片描述

绘制决策边界等价于绘制 $θ^Tx=0$ 线。

x = load('ex4x.dat');
y = load('ex4y.dat');
m = length(y);
mu = mean(x);
xx = x;
sigma = std(x);
x = (x - mean(x))./std(x);
x = [ones(m,1),x] ;
n = size(x,2);
% find返回满足指定条件的行的索引
pos = find ( y == 1 ) ; neg = find ( y == 0 ) ;
plot ( xx ( pos , 1 ) , xx ( pos , 2 ) , '+' ) ; 
hold on
plot ( xx ( neg , 1 ) , xx ( neg , 2 ) , ' o ' )
xlabel('Exam1 score');
ylabel('Exam2 score');
legend('Admitted','Not admitted');

MaxIter = 1500;
theta = zeros(size(x(1,:)))';
e = 1e-6;
H = zeros(n,n);
g = @(x) 1./(1+exp(-x));
for i = 1:MaxIter
    z = x * theta;
    h = g(z); %logistic model
    L_theta(i,1) = -(1/m)*sum(y.*log(h)+(1-y).*log(1-h)); %log likelihood function
    delta_L = (1/m)*x'*(h-y); %calculate gradient   
    %calculate Hessian matrix
    H = (1/m).*x' * diag(h) * diag(1-h) * x;
    %update L and theta
    if (i > 1) && (abs(L_theta(i,1) - L_theta(i-1,1)) <= e )
        break;
    end
    theta = theta - H^(-1)*delta_L;
    store(i,:) = [theta',L_theta(i,1)];
end
x_axis = x(:,2)*sigma(1) + mu(1);
y_axis = (-theta(1,1).*x(:,1) - theta(2,1).*x(:,2))/theta(3,1);
y_axis = y_axis*sigma(2) + mu(2);
plot(x_axis, y_axis,'-');
figure
plot(1:i-1,store,'-');
legend('\theta_0','\theta_1','\theta_2','L{(\theta)}');
xlabel('iter value');
ylabel('value');