Exercise: Logistic Regression and Newton's Method

題目地址：

Exercise: Logistic Regression

題目概要：某個高中有80名學生，其中40名得到了大學的錄取，40名沒有被錄取。x中包含80名學生兩門標准考試的成績，y中包含學生是否被錄取（1代表錄取、0代表未錄取）。

過程：

1、加載試驗數據，並為x輸入添加一個偏置項。

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

2、繪制數據分布

 % find returns the indices of the
% rows meeting the specified condition
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')

3、Newton's method

先回憶一下logistic regression的假設：

因為matlab中沒有sigmoid函數，因此我們用niline定義一個：

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

再來看看我們定義的cost function J(θ):

我們希望使用Newton's method來求出cost function J(θ)的最小值。回憶一下Newton's method中θ的迭代規則為：

在logistic regression中，梯度和Hessian的求法分別為：

需要注意的是，上述公式中的寫法為向量形式。

其中，為n+1 * 1的向量，為n+1 * n+1的矩陣。

和為標量。

實現

 按照上述Newton's method所述的方法逐步實現。代碼如下

function [theta, J ] = newton( x,y )
%NEWTON Summary of this function goes here
% Detailed explanation goes here
m = length(y);
theta = zeros(3,1);
g = inline('1.0 ./ (1.0 + exp(-z))'); 
pos = find(y == 1);
neg = find(y == 0);
J = zeros(10, 1); 
for num_iterations = 1:10
    %計算實際輸出
    h_theta_x = g(x * theta);
    %將y=0和y=1的情況分開計算后再相加，計算J函數
    pos_J_theta = -1 * log(h_theta_x(pos));
    neg_J_theta = -1 *log((1- h_theta_x(neg)));
    J(num_iterations) = sum([pos_J_theta;neg_J_theta])/m;
    %計算J導數及Hessian矩陣
    delta_J = sum(repmat((h_theta_x - y),1,3).*x);
    H = x'*(repmat((h_theta_x.*(1-h_theta_x)),1,3).*x);
    %更新θ
    theta = theta - inv(H)*delta_J';
end
% now plot J
% technically, the first J starts at the zero-eth iteration
% but Matlab/Octave doesn't have a zero index
figure;
plot(0:9, J(1:10), '-')
xlabel('Number of iterations')
ylabel('Cost J')
end

PS：在練習的附錄答案中給出的直接計算J函數的代碼很優雅：

J(i) =(1/m)*sum(-y.*log(h) - (1-y).*log(1-h));

在matlab中調用：

[theta J] = newton(x,y);

可以得到輸出：

θ為：

J函數輸出圖像為（令迭代次數為10的情況下）：

可以看到，其實在迭代至第四次時，就已經收斂了。

我們可以輸出一下admitted和unadmitted的分界線，代碼如下：

 plot(x(pos, 2), x(pos,3), '+'); hold on
plot(x(neg, 2), x(neg, 3), 'o')
ylabel('exam 2 scores')
xlabel('exam 1 scores')
plot(x(:,2), (theta(1) - x(:,2)*theta(2))/theta(3), '-');
legend('admitted', 'unadmitted','decision boundary');

效果如圖：