% Exercise 4 -- Logistic Regression
clear all; close all; clc
x = load(‘E:\workstation\data\ex4x.dat‘);
y = load(‘E:\workstation\data\ex4y.dat‘);
[m, n] = size(x);
% Add intercept term to x
x = [ones(m, 1), x];
% Plot the training data
% Use different markers for positives and negatives 分类分别计算每个的具体是允许还是不允许
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(‘Exam 1 score‘)
ylabel(‘Exam 2 score‘)
% Initialize fitting parameters
theta = zeros(n+1, 1);
% Define the sigmoid function
g = inline(‘1.0 ./ (1.0 + exp(-z))‘);
% Newton‘s method
MAX_ITR = 7;
J = zeros(MAX_ITR, 1);
for i = 1:MAX_ITR
% Calculate the hypothesis function
z = x * theta;
h = g(z);%转换成logistic函数
% Calculate gradient and hessian.
% The formulas below are equivalent to the summation formulas
% given in the lecture videos.
grad = (1/m).*x‘ * (h-y);%梯度的矢量表示法
H = (1/m).*x‘ * diag(h) * diag(1-h) * x;%hessian矩阵的矢量表示法
% Calculate J (for testing convergence)
J(i) =(1/m)*sum(-y.*log(h) - (1-y).*log(1-h));%损失函数的矢量表示法
theta = theta - H\grad;%是这样子的吗？
end
% Display theta
% Calculate the probability that a student with
% Score 20 on exam 1 and score 80 on exam 2
% will not be admitted
prob = 1 - g([1, 20, 80]*theta)
%画出分界面
% Plot Newton‘s method result
% Only need 2 points to define a line, so choose two endpoints
plot_x = [min(x(:,2))-2, max(x(:,2))+2];
% Calculate the decision boundary line，plot_y的计算公式见博客下面的评论。
plot_y = (-1./theta(3)).*(theta(2).*plot_x +theta(1));
plot(plot_x, plot_y)
legend(‘Admitted‘, ‘Not admitted‘, ‘Decision Boundary‘)
hold off
% Plot J
figure
plot(0:MAX_ITR-1, J, ‘o--‘, ‘MarkerFaceColor‘, ‘r‘, ‘MarkerSize‘, 8)
xlabel(‘Iteration‘); ylabel(‘J‘)
% Display J
结果：
logistic regression练习
这里给出的训练样本的特征为80个学生的两门功课的分数，样本值为对应的同学是否允许被上大学，如果是允许的话则用’1’表示，否则不允许就用’0’表示，这是一个典型的二分类问题。在此问题中，给出的80个样本中正负样本各占40个。而这节采用的是logistic regression来求解，该求解后的结果其实是一个概率值，当然通过与0.5比较就可以变成一个二分类问题了。
实验基础：
在logistic regression问题中，logistic函数表达式如下：
这样做的好处是可以把输出结果压缩到0~1之间。而在logistic回归问题中的损失函数与线性回归中的损失函数不同，这里定义的为：
如果采用牛顿法来求解回归方程中的参数，则参数的迭代公式为：
其中一阶导函数和hessian矩阵表达式如下：
当然了，在编程的时候为了避免使用for循环，而应该直接使用这些公式的矢量表达式（具体的见程序内容）。
一些matlab函数：
find:
是找到的一个向量，其结果是find函数括号值为真时的值的下标编号。
inline:
构造一个内嵌的函数，很类似于我们在草稿纸上写的数学推导公式一样。参数一般用单引号弄起来，里面就是函数的表达式，如果有多个参数，则后面用单引号隔开一一说明。比如：g = inline(‘sin(alpha*x)‘,‘x‘,‘alpha‘)，则该二元函数是g(x,alpha) = sin(alpha*x)。
Min max 数组中最大最小元素
具体步骤：
1：加载数据和画图：
x = load(‘E:\workstation\data\ex4x.dat‘);
y = load(‘E:\workstation\data\ex4y.dat‘);
[m, n] = size(x);计算数据行列
% Add intercept term to x
x = [ones(m, 1), x]; 将第一列数据变为1
% Plot the training data
% Use different markers for positives andfigure
pos = find(y); neg = find(y ==0)%找到对应允许和不允许的分别画标记；
plot(x(pos, 2), x(pos,3), ‘+‘)
hold on
plot(x(neg, 2), x(neg, 3), ‘o‘)
hold on
xlabel(‘Exam 1 score‘)
ylabel(‘Exam 2 score‘)
Newton‘s Method
Recall that in logistic regression, the hypothesis function is
In our example, the hypothesis is interpreted as the probability that a driver will be accident-free, given the values of the features in x.
Matlab/Octave does not have a library function for the sigmoid, so you will have to define it yourself. The easiest way to do this is through an inline expression:
g = inline(‘1.0 ./ (1.0 + exp(-z))‘);
% Usage: To find the value of the sigmoid
% evaluated at 2, call g(2)
The cost function is defined as
Our goal is to use Newton‘s method to minimize this function. Recall that the update rule for Newton‘s method is
In logistic regression, the gradient and the Hessian are
Note that the formulas presented above are the vectorized versions. Specifically, this means that , , while and are scalars.
Implementation
Now, implement Newton‘s Method in your program, starting with the initial value of . To determine how many iterations to use, calculate for each iteration and plot your results as you did in Exercise 2. As mentioned in the lecture videos, Newton‘s method often converges in 5-15 iterations. If you find yourself using far more iterations, you should check for errors in your implementation.
After convergence, use your values of theta to find the decision boundary in the classification problem. The decision boundary is defined as the line where
which corresponds to
% Initialize fitting parameters
theta = zeros(n+1, 1); 初始化θ值为0。
% Define the sigmoid function
g = inline(‘1.0 ./ (1.0 + exp(-z))‘); 定义归一化函数使g在【0,1】之间inline为内联函数
% Newton‘s method
MAX_ITR = 7;
J = zeros(MAX_ITR, 1);牛顿法具体做法
for i = 1:MAX_ITR
% Calculate the hypothesis function
z = x * theta; % 有公式定义值 h = g(z);
h = g(z);
grad = (1/m).*x‘ * (h-y); 梯度下降
H = (1/m).*x‘ * diag(h) * diag(1-h) * x;%牛顿环矢量表示
J(i) =(1/m)*sum(-y.*log(h) - (1-y).*log(1-h));随时函数矢量表示
theta = theta - H\grad; 求出θ这样可以画出直线y=θ0+θ1x曲线从而判断具体可以被允许
end
% Display theta
3：下面找到分界面尽量将一样的放在一方然后之后的可以根据分界面预测是否被允许、
% Calculate the probability that a student with
% Score 20 on exam 1 and score 80 on exam 2
% will not be admitted
prob = 1 - g([1, 20, 80]*theta) %测试一下数据为【1,20,80】这个人
% Plot Newton‘s method result
% Only need 2 points to define a line, so choose two endpoints
plot_x = [min(x(:,2))-2, max(x(:,2))+2]; 第二列即学期最大最小值 +-2是原来的线更加延伸，
直接令logistic回归的值为0.5，则可以得到e的指数为0，即：
theta(1)*1+theta(2)*plot_x+theta(3)*plot_y=0,解出plot_y即可。
plot_y = (-1./theta(3)).*(theta(2).*plot_x +theta(1));
plot(plot_x, plot_y)
legend(‘Admitted‘, ‘Not admitted‘, ‘Decision Boundary‘)
hold off
% Plot J 损失函数值和迭代次数之间的曲线：
figure
plot(0:MAX_ITR-1, J, ‘o--‘, ‘MarkerFaceColor‘, ‘r‘, ‘MarkerSize‘, 8)
xlabel(‘Iteration‘); ylabel(‘J‘)
% Display J