逻辑斯蒂（logistic）回归深入理解、阐述与实现

逻辑斯蒂（logistic）回归深⼊理解、阐述与实现
第⼀节中说了，logistic 回归和线性回归的区别是：线性回归是根据样本X各个维度的Xi的线性叠加（线性叠加的权重系数wi就是模型的参数）来得到预测值的Y，然后最⼩化所有的样本预测值Y与真实值y'的误差来求得模型参数。

我们看到这⾥的模型的值Y是样本X各个维度的Xi的线性叠加，是线性的。

Y=WX (假设W>0),Y的⼤⼩是随着X各个维度的叠加和的⼤⼩线性增加的，如图（x为了⽅便取1维）：
然后再来看看我们这⾥的logistic 回归模型，模型公式是：，这⾥假设W>0,Y与X各维度叠加和（这⾥都是线性叠加W）的图形关系，如图（x为了⽅便取1维）：
我们看到Y的值⼤⼩不是随X叠加和的⼤⼩线性的变化了，⽽是⼀种平滑的变化，这种变化在x的叠加和为0附近的时候变化的很快，⽽在很⼤很⼤或很⼩很⼩的时候，X叠加和再⼤或再⼩，Y值的变化⼏乎就已经很⼩了。

当X各维度叠加和取⽆穷⼤的时候，Y趋近于1，当X各维度叠加和取⽆穷⼩的时候，Y趋近于0.
这种变量与因变量的变化形式就叫做logistic变化。

（注意不是说X各个维度和为⽆穷⼤的时候，Y值就趋近1，这是在基于W>0的基础上，（如果W<0,n那么Y趋近于0）⽽W是根据样本训练出来，可能是⼤于0，也可能是⼩0，还可能W1>0，W2<0…所以这个w值是样本⾃动训练出来的，也因此不是说你只要x1，x2，x3…各个维度都很⼤，那么Y值就趋近于1，这是错误的。

凭直觉想⼀下也不对，因为你连样本都还没训练，你的模型就有⼀个特点：X很⼤的时候Y就很⼤。

这种强假设肯定是不对的。

因为可能样本的特点是X很⼤的时候Y就很⼩。

）
所以我们看到，在logistic回归中，X各维度叠加和（或X各维度）与Y不是线性关系，⽽是logistic关系。

⽽在线性回归中，X各维度叠加和就是Y，也就是Y与X就是线性的了。

X各维度叠加和Y的关系不只是这⼀种，还可能是其他的⽐如：
为什么变量与因变量要选⽤logistic关系呢，因为这⾥（1）我们需要Y代表的是概率即Y∈（0,1）。

（2）我们需要X各维度叠加和在0附近变化幅度⽐较⼤，并且是⾮线性的变化。

⽽在很⼤或很⼩的时候，⼏乎不变化，这是基于概率的⼀种认识与需要。

感性的⼀个例⼦，想想你学习努⼒的程度与从60分提⾼到80分和80提⾼到100分并不是线性的。

（3）这个关系的公式要在之后形成的cost function是凸函数。

所以就选择了logistic。

前⾯已经说了，我们使⽤logistic回归是⽤于⼆分类问题（y只有两个值A,B，也可以写成1和0，这都没关系），回归模型得到的结果不是预测样本X对应的y值（注意下，在logistic回归这⾥我们⼩写y表⽰某个样本Xi的类别，⽽⼤写Y或Y(Xi)表⽰logistic回归模型对某个样本Xi预测为1的概率。

其实这⾥最好把Y⽤其他字母表⽰，以免混淆，但是已经这⾥写了，以后注意。

），⽽是y=1的概率或y=0的概率。

我们假设y=1的概率公式是：，那么y=0的概率就是。

（注意我们也可以y=0的概率公式为前⾯那⼀个，这⾥是任意的。

这⾥不同的结果只是最终的W参数不同罢了。

因为我们最终的W是训练出来的，不管怎么样，模型都会表现出样本的特点来。

只是我们习惯了把Y(X)当成y=1的logistic模型映射的概率）
还要注意这⾥我们不是对⼀个Xi都要分别预测出来y=1的概率和y=0的概率。

⽽是对于⼀个Xi，如果它的yi=1，那么我们就⽤这个公式映射所对应的概率，如果对于⼀个Xi，如果它的yi=0，那么我们就⽤这个公式映射所对应的概率。

都是根据yi的值映射出来⼀个概率。

因为我们的Y是概率，我们不能利⽤最⼩误差等，我们这⾥⽤的是极⼤化所有样本的对数似然函数：。

yi表⽰的是Xi真实所属的类别（1或0）。

L（W）就是cost function。

这⾥的⽬标函数值和W有关系，也就是X各维度线性叠加的那些权重有关系。

那么这些权重的物理意义是什么呢？就是X各维度与Y值关系的那个⽅向，听起来可能很抽象，现在看⼀下具体例⼦（X为⼆维,但考虑将常数项变成齐次后X就是三维了，第⼀维是1，因此W也是三维的，第⼀维是常数项，⼤⼩不影响那个⽅向，主要考虑后⾯的两个值的关系），⽐如：
当W=[-15 0.2 0.2]，图形为：（⼀种45°⽅向的关系）
当W=[-15 -0.2 -0.2]，图形为：（⼀种-45°⽅向的关系）
当W=[-15 0.2 0.01]，图形为：（⼀种0°⽅向（沿x轴形成的logistic关系）的关系）
当W=[-15 0.01 0.2]，图形为：（⼀种90°⽅向（沿y轴形成的logistic关系）的关系）
下⾯我们对L(W)求极值。

L（W）是negtive 的似然函数。

只有是negtive的，才能⽤梯度下降法求最⼩值，如果不是negtive的，就要⽤梯度上升法求最⼤值了，我们⼀般不⽤那个。

还要注意我们的代价函数不管是最⼩⼆乘法、误差均⽅、似然函数等，都要求平均，就是前⾯加上1/m.利⽤梯度下降法求得的迭代公式是：，其中wj代表的是模型参数向量W的第j个元素。

α代表的是学习速率，yi表⽰的是第i个样本向量的真实标签（label），也就是第i个样本向量所属的类别（0或1），Y（Xi）表⽰的是回归模型预测的第i个样本向量为1的概率。

xij表⽰的第i个样本向量Xi的第j个元素。

⼩⼼不要忘了Σ（i=1：m）。

还要注意因为梯度下降法梯度是⽬标函数分别对模型参数向量W的第每⼀个元素求偏导，所以这⾥W的第每⼀个元素的值是⼀个⼀个求出来的。

当然在matlab中，虽然是按向量导⼊计算的数据，但是本质上还是⼀个⼀个计算的w的每个值，⽽不是直接求的这个向量是多少。

后来⼜看了看，前⾯说的这就句话不对，虽然求偏导是分别对W的每⼀维度求，但是这个L（w）每⼀次迭代的梯度向量是可以⼀次求出来的，也是可以在公式⾥直接表达出来
我们注意到这个迭代公式和线性回归的迭代公式形式是⼀样的，只是在线性回归中Y(Xi)为Xi的各维度线性叠加即WXi=WXi，⽽这⾥Xi的各维度线性叠加WXi以后,还要进⾏⼀次⾮线性映射（即logistic映射），⾮线性的映射到（0,1）之间Y(Xi)。

所以也可以认为logstic回归也是处理的线性问题，即也是求的各维度线性叠加的权重系数，只是求得了各维度线性叠加和以后，不是与Xi所属的类别进⾏⽐较，⽽是⾮线性映射到所属类别的概率，然后最⼤似然函数法求模型参数（也就是那个各维度线性叠加的权重，⽽后⾯的⾮线性映射logstic那⾥没有参数）。

下⾯说⼀个例⼦：训练样本的特征为80个学⽣的两门功课的分数，样本值yi为对应的同学是否允许被上⼤学。

训练好以后，⽤训练好的模型根据某个学⽣的成绩预测是否被允许上⼤学？数据中给出被允许上⼤学为1，不被允许上⼤学为0。

程序代码（梯度下降法）：
1: clear all; close all; clc
2: x = load('ex4x.dat'); %每⼀⾏是⼀个样本
3: y = load('ex4y.dat');
4: [m, n] = size(x);
5: sample_num = m;
6: x = [ones(m, 1), x]; %x增加⼀维。

因为前⾯⽂字或前⼏节说了
7: % Plot the training data
8: % Use different markers for positives and negatives
9: figure;
10: pos = find(y == 1); neg = find(y == 0);%pos 和neg 分别是 y元素=1和0的所在的位置序号组成的向量
11: plot(x(pos, 2), x(pos,3), '+')%⽤+表⽰那些yi=1所对应的样本
12: hold on
13: plot(x(neg, 2), x(neg, 3), 'o')
14: hold on
15: xlabel('Exam 1 score')
16: ylabel('Exam 2 score')
17: itera_num=500;%迭代次数
18: g = inline('1.0 ./ (1.0 + exp(-z))'); %这个就相当于制造了⼀个function g（z）=1.0 ./ (1.0 + exp(-z))
19: plotstyle = {'b', 'r', 'g', 'k', 'b--', 'r--'};
20: figure;%建⽴新的窗⼝
21: alpha = [ 0.0009, 0.001,0.0011,0.0012,0.0013 ,0.0014 ];%下⾯就分别⽤这⼏个学习速率看看哪个更好
22:for alpha_i = 1:length(alpha) %alpha_i是1,2，...6，表⽰的是学习速率向量和曲线格式向量的坐标：alpha(alpha_i)，plotstyle(alpha_i)
23: theta = zeros(n+1, 1);%thera表⽰样本Xi各个元素叠加的权重系数，这⾥以向量形式表⽰，且初始都为0，三维向量
24: J = zeros(itera_num, 1);%J是个100*1的向量，第n个元素代表第n次迭代cost function的值（下⾯⽤negtive 的对数似然函数，
25: %因为是negtive 的，所以是求得极⼩值）
26:for i = 1:itera_num %计算出某个学习速率alpha下迭代itera_num次数后的参数
27: z = x * theta;%这个z是⼀个列向量，每⼀个元素是每⼀个样本Xi的线性叠加和，因为X是所有的样本，因此这⾥不是⼀个⼀个样本算的，
28: %⽽是所有样本⼀块算的，因此z是⼀个包含所有样本Xi的线性叠加和的向量。

在公式中，是单个样本表⽰法，⽽在matlab中都是所有的样本⼀块来。

29: h = g(z);%这个h就是样本Xi所对应的yi=1时，映射的概率。

如果⼀个样本Xi所对应的yi=0时，对应的映射概率写成1-h。

30: J(i) =(1/sample_num).*sum(-y.*log(h) - (1-y).*log(1-h));%损失函数的⽮量表⽰法这⾥Jtheta是个100*1的列向量。

31: grad = (1/sample_num).*x'*(h-y);%在程序中X是包括所有的X样本，公式中只能是每⼀个样本Xi，所以不能直接把公式照搬，所以要认真看看，代码中相应改变⼀下。

33: theta = theta - alpha(alpha_i).*grad;
34: end
35: plot(0:itera_num-1, J(1:itera_num),char(plotstyle(alpha_i)),'LineWidth', 2)
%此处⼀定要通过char函数来转换因为包⽤（）索引后得到的还是包cell，
36: %所以才要⽤char函数转换，也可以⽤{}索引，这样就不⽤转换了。

37: %⼀个学习速率对应的图像画出来以后再画出下⼀个学习速率对应的图像。

38: hold on
39: if(1 == alpha(alpha_i)) %通过实验发现alpha为0.0013 时效果最好，则此时的迭代后的theta值为所求的值
40: theta_best = theta;
41: end
42: end
43: legend('0.0009', '0.001','0.0011','0.0012','0.0013' ,'0.0014');%给每⼀个线段格式标注上
44: xlabel('Number of iterations')
45: ylabel('Cost function')
46: prob = g([1, 20, 80]*theta);
%把[1, 20, 80]*theta这个带⼊g（z）得到的就是在exam1=20、exam2=80的条件下，通过的概率（也就是Y=1）的概率是多少。

47: %画出分界⾯
48: % Only need 2 points to define a line, so choose two endpoints
49: plot_x = [min(x(:,2))-2, max(x(:,2))+2];%两个点就可以画出来直线，这⾥这么取x1坐标是让直接的视野更容易包含那些样本点。

50: plot_y = (-1./theta(3)).*(theta(2).*plot_x +theta(1));%分界⾯怎么画呢？问题也就是在x1，x2坐标图中找到
%那些将x1，x2带⼊1/(1+exp(-wx))后，使其值>0.5的（x1,x2）
51: %坐标形成的区域，因为我们知道1/(1+exp(-wx))>0.5
52: %意味着该区域的（x1，x2）表⽰的成绩允许上⼤学的概率>0.5,那么其他的区域就是不被允许上⼤学，那么1/(1+exp(-wx))=0.5解出来的⼀个关
%于x1，x2的⽅程就是那个分界⾯。

53: %我们解出来以后发现，这个⽅程是⼀个直线⽅程：w（2）x1+w（3）x2+w（1）=0
54: %注意我们不能因为这个分界⾯是直线，就认为logistic回归是⼀个线性分类器，注意logistic回归不是⼀个分类器，他没有分类的功能，
%这个logistic回归是⽤来
55: %预测概率的，这⾥具有分类功能是因为我们硬性规定了⼀个分类标准：把>0.5的归为⼀类，<0.5的归于另⼀类。

这是⼀个很强的假设，
%因为本来我们可能预测了⼀个样本
56: %所属某个类别的概率是0.6，这是⼀个不怎么⾼的概率，但是我们还是把它预测为这个类别，只因为它>0.5.所以最后可能logistic回归加上这
%个假设以后形成的分类器
57: %的分界⾯对样本分类效果不是很好，这不能怪logistic回归，因为logistic回归本质不是⽤来分类的，⽽是求的概率。

58: figure;
59: plot(x(pos, 2), x(pos,3), '+')%把分界⾯呈现在样本点上⾯，所以还要把样本点画出来。

60: hold on
61: plot(x(neg, 2), x(neg, 3), 'o')
62: hold on
63: plot(plot_x, plot_y)
64: legend('Admitted', 'Not admitted', 'Decision Boundary')
65: hold off
66:
当学习速率取不同的值时，迭代次数与cost function形成的不同的曲线如图所⽰：
当学习速率为0.0014时候，我们看到图像开始振荡，说明学习速率过⼤了。

当我们将logistic回归应⽤到分类器中时，如果我们假设p（Y(x)|x）>0.5归为⼀类，p（Y(x)|x）<0.5归为另⼀类。

那么分类器分界⾯如图：当我们看到这个图时，我们是不是有种感觉：这个logistic回归⽐线性回归复杂，为什么结果得到的那么差？
先说⼀句，我们不要⽤这个分界⾯来评价logistic回归模型的好坏这个⼀会再说，我们先说这个分界⾯公式是怎么产⽣的。

分界⾯怎么画呢？问题也就是在x1，x2坐标图中找到那些将x1，x2带⼊1/(1+exp(-wx))后，使其值>0.5的（x1,x2）坐标形成的区域，因为我们知道1/(1+exp(-wx))>0.5意味着该区域的（x1，x2）表⽰的成绩允许上⼤学的概率>0.5,那么其他的区域就是不被允许上⼤学，那么
1/(1+exp(-wx))=0.5解出来的⼀个关于x1，x2的⽅程就是那个分界⾯。

我们解出来以后发现，这个⽅程是⼀个直线⽅程：
w（2）x1+w（3）x2+w（1）=0 注意我们不能因为这个分界⾯是直线，就认为logistic回归是⼀个线性分类器，注意logistic回归不是⼀个分类器，他没有分类的功能，这个logistic回归是⽤来预测概率的，这⾥具有分类功能是因为我们硬性规定了⼀个分类标准：把>0.5的归为⼀类，<0.5的归于另⼀类。

这是⼀个很强的假设，因为本来我们可能预测了⼀个样本所属某个类别的概率是0.6，这是⼀个不怎么⾼的概率，但是我们还是把它预测为这个类别，只因为它>0.5.所以最后可能logistic回归加上这个假设以后形成的分类器
的分界⾯对样本分类效果不是很好，这不能怪logistic回归，因为logistic回归本质不是⽤来分类的，⽽是求的概率。

⽜顿法
注意⼀下，在这个例⼦中，我们要求G，我们不是利⽤这个公式，这个太⿇烦了，⽽是直接在刚才⽬标函数⼀次求导的基础上继续求导，⼆次导为：
matlab程序中的形式具体看程序。

⽜顿法程序clear all; close all; clx = load('ex4x.dat');
y = load('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);%是可以以向量形式给g（）函数输⼊的，不是因为所有函数都可以这么输，⽽是这个函数⾥⾯有个exp（）函数，
%是可以的。

z的每⼀维在exp（）中分别计算，对应h⾥的每⼀维。

% 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;%因为在matlab中，是以所有样本⼀块作为X输⼊的，所以这⾥的h是⼀个
%向量，⾥⾯的元素是每个样本映射的logistic回归模型的概率，这样两个对⾓矩阵乘出来，还是⼀个对⾓矩阵，每个对⾓
%元是每个样本的hi*（1-hi）。

这样做就是适应公式在matlab下的表现形式，是正确的，我验算了。

% Calculate J (for testing convergence)
J(i) =(1/m)*sum(-y.*log(h) - (1-y).*log(1-h));
theta = theta - H\grad;%是这样的，matlab⾥ \ 是左除，/ 是右除。

如果是数字计算，则左除和右除是等效的，例如 3/2 = 2\3。

%⽽对于矩阵运算，则⼆者不等效。

矩阵除法在 matlab ⾥定义为矩阵求逆后相乘。

例如 A的逆矩阵是 A1，则 B/A = B*A1，A\B = A1*B。

%矩阵乘法不满⾜交换律，因此需要有左右除法之分。

矩阵求逆的命令是 inv矩阵乘法不满⾜交换律，因此需要有左右除法之分。

end
% Display theta
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
J
我们可以看到利⽤⽜顿法是⽐梯度下降法收敛快。

这个解释的原因很多，在这⾥就不解释了。

完。