文章目录
一、逻辑回归算法简介
逻辑回归(Logistic regression)是最经典的分类算法之一,一般入门机器学习也是从逻辑回归和线性回归这两个算法入手,线性回归的推导,已经在梯度下降法哪里涉及了,本文主要是关于逻辑回归的推导以及使用梯度下降法和牛顿法来求解逻辑回归。
二、逻辑回归的数学思想
1.构造预测函数
前面说了逻辑回归是用于二分类问题,(即输出是有两种1或者0)所以逻辑回归的输出值应该也只有两个值(这两个值分别代表两个类别),那么怎么才能找到一个函数的输出值只有两个值呢?
引入一个新的函数概念:Sigmoid函数,这个函数的定义如下:
g
(
x
)
=
1
1
+
e
−
x
g(x) = \frac{1}{{1 + {e^{ - x}}}}
g(x)=1+e−x1
这个函数的值域应该是0到1,自变量是
x
x
x,所以可以绘制此函数的图像:
可以这样理解,当x的值大于0,取
g
(
x
)
g(x)
g(x)的值为1,当
x
x
x的值小于0,取
g
(
x
)
g(x)
g(x)的值为0,这样就可以达到输出为0或者1的目的。
在找到了分类所用函数之后,假设需要解决的问题是将下面的数据进行分类:
自变量的个数有两个:
x
1
、
x
2
x_{1}、x_{2}
x1、x2,需要找到一条直线,能够合理的将上述两类点进行分开,即
h
(
x
)
h(x)
h(x),
h
(
x
)
h(x)
h(x)表示这条直线。我们的任务就是找到这条直线的参数
θ
\theta
θ,结合这条直线和Sigmoid函数得到:
h
θ
(
x
)
=
g
(
θ
T
x
)
=
1
1
+
e
−
θ
T
x
{h_\theta }(x) = g({\theta ^T}x) = \frac{1}{{1 + {e^{ - {\theta ^T}x}}}}
hθ(x)=g(θTx)=1+e−θTx1
上述就是逻辑回归的预测函数,
h
(
x
)
h(x)
h(x)的值应该在0到1之间,它是结果取1的概率值。
综上所述我们便得到了对于输入数据进行分类结果的类别1和类别2 的概率分别为
P
(
y
=
1
∣
x
;
θ
)
=
h
θ
(
x
)
P(\left. {y = 1} \right|x;\theta ) = {h_\theta }(x)
P(y=1∣x;θ)=hθ(x)
P
(
y
=
1
∣
x
;
θ
)
=
1
−
h
θ
(
x
)
P(\left. {y = 1} \right|x;\theta ) = 1 - {h_\theta }(x)
P(y=1∣x;θ)=1−hθ(x)
2.构造损失函数
在得到逻辑回归的数学形式之后,剩下的就是如何去求解模型中的参数
θ
\theta
θ,关于如何求解这个参数,在机器学习中,我们常常会引入某个损失函数,使得整个损失函数的值最小的参数就是我们需要求解的参数。由于逻辑回归是一个概率模型,我们可以使用极大似然估计法来求解模型的参数。这里先来看一下似然估计函数为:
L
(
w
)
=
∏
[
h
(
x
i
)
]
y
i
[
1
−
h
(
x
i
)
]
1
−
y
i
L(w) = \prod {{{[h({x_i})]}^{{y_i}}}{{[1 - h({x_i})]}^{1 - {y_i}}}}
L(w)=∏[h(xi)]yi[1−h(xi)]1−yi
对上述的式子两边都取对数转换,得到对数似然函数:
ln
(
L
(
θ
)
)
=
∑
[
y
i
ln
h
(
x
i
)
+
(
1
−
y
i
)
ln
(
1
−
h
(
x
i
)
)
]
\ln (L(\theta )) = \sum {[{y_i}\ln h({x_i}) + (1 - {y_i})\ln (1 - h({x_i}))]}
ln(L(θ))=∑[yilnh(xi)+(1−yi)ln(1−h(xi))]
对于上述式子的右边我们将多项式合并可以得到:
ln
(
L
(
θ
)
)
=
∑
[
y
i
ln
h
(
x
i
)
1
−
h
(
x
i
)
+
ln
(
1
−
h
(
x
i
)
)
]
\ln (L(\theta )) = \sum {[{y_i}\ln \frac{{h({x_i})}}{{1 - h({x_i})}} + \ln (1 - h({x_i}))]}
ln(L(θ))=∑[yiln1−h(xi)h(xi)+ln(1−h(xi))]
到这里需要对
h
(
x
i
)
1
−
h
(
x
i
)
{\frac{{h({x_i})}}{{1 - h({x_i})}}}
1−h(xi)h(xi)进行一个变化,仔细看,
h
(
x
i
)
h(x_{i})
h(xi)表示函数取值类别为1的概率,即有:
h
(
x
i
)
=
1
1
+
e
−
θ
T
x
h({x_i}) = \frac{1}{{1 + {e^{ - {\theta ^T}x}}}}
h(xi)=1+e−θTx1
得到:
θ
T
x
i
=
ln
h
(
x
i
)
1
−
h
(
x
i
)
{\theta ^T}{x_i} = \ln \frac{{h({x_i})}}{{1 - h({x_i})}}
θTxi=ln1−h(xi)h(xi)
将
θ
T
x
\theta^{T}x
θTx的表达式带入对数似然函数中可以得到:
ln
(
L
(
θ
)
)
=
∑
[
y
i
ln
(
θ
T
x
)
+
ln
(
1
−
h
(
x
i
)
)
]
\ln (L(\theta )) = \sum {[{y_i}\ln (\theta^{T}x)+ \ln (1 - h({x_i}))]}
ln(L(θ))=∑[yiln(θTx)+ln(1−h(xi))]
而对于右边的
l
n
(
1
−
h
(
x
i
)
)
ln(1-h(x_{i}))
ln(1−h(xi))我们也是可以推导出来的。
即
ln
(
1
−
h
(
x
i
)
)
=
ln
(
1
−
1
1
+
e
−
θ
T
x
)
=
ln
(
1
e
θ
T
x
+
1
)
\ln(1-h(x_{i}))=\ln(1-\frac{1}{1+e^{-\theta^{T}x}})=\ln(\frac{1}{e^{\theta^{T}x}+1})
ln(1−h(xi))=ln(1−1+e−θTx1)=ln(eθTx+11)
最终可以得到对数似然函数为:
ln
(
L
(
θ
)
)
=
∑
[
y
i
ln
(
θ
T
x
)
−
ln
(
e
θ
T
x
+
1
)
]
\ln(L(\theta))=\sum[y_{i}\ln(\theta^{T}x)-\ln(e^{\theta^{T}x}+1)]
ln(L(θ))=∑[yiln(θTx)−ln(eθTx+1)]
我们上面说过在机器学习中常会用损失函数最小化来求解参数,在逻辑回归中取整个数据集上的平均对数损失函数作为其学习参数的损失函数,即:
J
(
θ
)
=
−
1
N
ln
(
L
(
θ
)
)
J(\theta)=-\frac{1}{N}\ln(L(\theta))
J(θ)=−N1ln(L(θ))
于是我们就可以把原本的极大化似然函数转换成最小化损失函数
J
(
θ
)
J(\theta)
J(θ)。
3.求解参数
求解参数主要任务就是求解逻辑回归中的
θ
\theta
θ,而且这里的
θ
\theta
θ还应该是三个值,因为:
θ
T
x
=
θ
1
x
1
+
θ
2
x
2
+
θ
\theta^{T} x=\theta_{1}x_{1}+\theta_{2}x_{2}+\theta
θTx=θ1x1+θ2x2+θ
这里的
θ
1
\theta_{1}
θ1和
θ
2
\theta_{2}
θ2表示自变量
x
1
x_{1}
x1和
x
2
x_2
x2的系数,
θ
\theta
θ表示常数项。
求解参数就是求解这三个参数的值。由于损失函数是光滑的凸函数,因此多种方法都可以使用,常用的就是梯度下降法和牛顿法了。
3.1.梯度下降法优化
首先介绍梯度下降法求解参数,在梯度下降法中,需要计算损失函数的一阶导函数,得到损失函数
J
(
θ
)
J(\theta)
J(θ)的一阶导函数为:
∂
J
(
θ
)
∂
θ
=
1
N
∑
(
h
(
x
i
)
−
y
i
)
x
i
\frac{{\partial J(\theta )}}{{\partial \theta }} = \frac{1}{N}\sum {(h({x_i}) - {y_i}){x_i}}
∂θ∂J(θ)=N1∑(h(xi)−yi)xi 更新参数的迭代式为:
θ
i
k
+
1
=
θ
i
k
−
α
∂
J
(
θ
)
∂
θ
i
\theta ^{k+1}_i=\theta^{k}_{i}-\alpha \frac{{\partial J(\theta )}}{{\partial {\theta _i}}}
θik+1=θik−α∂θi∂J(θ)
其中
k
k
k为迭代次数。这些都是梯度下降法的基本原理,详细文档参考:
机器学习笔记-梯度下降法
最终可以设置终止条件为迭代次数或者绝对误差小于阈值。
3.2.牛顿法优化
牛顿法与梯度下降法相比,是损失函数的二阶迭代,对损失函数作二阶泰勒展开,然后再求损失函数的二阶偏导数,对损失函数的二阶偏导数,其实就是对下面的式子再进行求导:
∂
J
(
θ
)
∂
θ
=
1
N
∑
(
h
(
x
i
)
−
y
i
)
x
i
\frac{{\partial J(\theta )}}{{\partial \theta }} = \frac{1}{N}\sum {(h({x_i}) - {y_i}){x_i}}
∂θ∂J(θ)=N1∑(h(xi)−yi)xi
∑
(
h
(
x
i
)
−
y
i
)
x
i
=
∑
(
1
1
+
e
−
θ
T
x
−
y
i
)
x
i
\sum(h(x_i)-y_i)x_i=\sum(\frac{1}{1+e^{-\theta^{T}x}}-y_i)x_i
∑(h(xi)−yi)xi=∑(1+e−θTx1−yi)xi
设
g
(
x
)
=
∑
(
e
θ
T
x
i
1
+
e
θ
T
x
i
−
y
i
)
x
i
g(x)=\sum(\frac{e^{\theta^{T}x_i}}{1+e^{\theta^{T}x_i}}-y_i)x_i
g(x)=∑(1+eθTxieθTxi−yi)xi,则
g
′
(
x
)
g'(x)
g′(x)为:
g
′
(
x
)
=
∑
(
x
i
e
θ
T
x
i
(
1
+
e
θ
T
x
)
−
x
i
e
θ
T
x
i
e
θ
T
x
i
(
1
+
e
θ
T
x
i
)
2
)
x
i
=
∑
x
i
(
e
θ
T
x
i
1
+
e
θ
T
x
i
)
x
i
=
∑
x
i
(
h
(
x
i
)
(
1
−
h
(
x
i
)
)
)
x
i
\begin{aligned} g'(x)&=\sum (\frac{x_ie^{\theta^{T}x_i}(1+e^{\theta^{T}x})-x_ie^{\theta^{T}x_i}e^{\theta^{T}x_i}}{(1+e^{\theta^{T}x_i})^{2}})x_i\\ &=\sum x_i(\frac{e^{\theta^{T}x_i}}{1+e^{\theta^{T}x_i}})x_i\\ &=\sum x_i(h(x_i)(1-h(x_i)))x_i \end{aligned}
g′(x)=∑((1+eθTxi)2xieθTxi(1+eθTx)−xieθTxieθTxi)xi=∑xi(1+eθTxieθTxi)xi=∑xi(h(xi)(1−h(xi)))xi
于是得到二阶偏导数为:
H
=
∂
2
J
(
θ
)
∂
θ
2
=
1
N
∑
x
i
(
h
(
x
i
)
(
1
−
h
(
x
i
)
)
)
x
i
H = \frac{{{\partial ^2}J(\theta )}}{{\partial {\theta ^2}}}=\frac{1}{N}\sum x_i(h(x_i)(1-h(x_i)))x_i
H=∂θ2∂2J(θ)=N1∑xi(h(xi)(1−h(xi)))xi
在得到黑塞矩阵后,便能够计算迭代方程式了
θ
k
+
1
=
θ
k
−
H
k
−
1
∂
J
k
(
θ
)
∂
θ
\theta^{k+1}=\theta^{k}-H^{-1}_k\frac{{\partial J_k(\theta )}}{{\partial \theta }}
θk+1=θk−Hk−1∂θ∂Jk(θ)
上面的
∂
J
k
(
θ
)
∂
θ
\frac{{\partial J_k(\theta )}}{{\partial \theta }}
∂θ∂Jk(θ)表示第k代的梯度。在得到两者的迭代表达式之后,我们便可以用代码来实现逻辑回归,由于Matlab中有求导的函数,所以我们采用Matlab来实现两种逻辑回归。
三、逻辑回归的应用
在这一节中,根据一个具体的题目来实现逻辑回归,先看训练的数据集:
从现有的数据集中看出,自变量一共有两个,我们需要找到一条合适的边界,能够完美的分离出这两类数据。
1.梯度下降法实现逻辑回归
代码:
主函数:
clc,clear % 清除命令
% 导入数据
tic
data = [ 34.623659624517 78.0246928153624 0
30.2867107682261 43.894997524001 0
35.8474087699387 72.9021980270836 0
60.1825993862098 86.3085520954683 1
79.0327360507101 75.3443764369103 1
45.0832774766834 56.3163717815305 0
61.1066645368477 96.5114258848962 1
75.0247455673889 46.5540135411654 1
76.0987867022626 87.420569719268 1
84.4328199612004 43.5333933107211 1
95.8615550709357 38.2252780579509 0
75.0136583895825 30.6032632342801 0
82.3070533739948 76.481963302356 1
69.3645887597094 97.7186919618861 1
39.5383391436722 76.0368108511588 0
53.9710521485623 89.207350137502 1
69.0701440628303 52.7404697301677 1
67.9468554771162 46.6785741067313 0
70.6615095549944 92.9271378936483 1
76.978783727475 47.5759636497553 1
67.3720275457088 42.8384383202918 0
89.6767757507208 65.7993659274524 1
50.534788289883 48.855811527642 0
34.2120609778679 44.2095285986629 0
77.9240914545704 68.9723599933059 1
62.2710136700463 69.9544579544759 1
80.1901807509566 44.8216289321835 1
93.114388797442 38.8006703371321 0
61.830206023126 50.2561078924462 0
38.7858037967942 64.9956809553958 0
61.379289447425 72.807887313171 1
85.4045193941164 57.0519839762712 1
52.1079797319398 63.1276237688172 0
52.0454047683183 69.4328601204522 1
40.2368937354511 71.1677480218488 0
54.6351055542482 52.2138858806112 0
33.9155001090689 98.8694357422061 0
64.1769888749449 80.9080605867082 1
74.7892529594154 41.5734152282443 0
34.1836400264419 75.2377203360134 0
83.9023936624915 56.3080462160533 1
51.5477202690618 46.8562902634998 0
94.4433677691785 65.5689216055905 1
82.3687537571392 40.6182551597062 0
51.0477517712887 45.82270145776 0
62.2226757612019 52.0609919483668 0
77.1930349260136 70.4582000018096 1
97.7715992800023 86.7278223300282 1
62.0730637966765 96.7688241241398 1
91.5649744980744 88.696292545466 1
79.9448179406693 74.1631193504376 1
99.2725269292572 60.9990309984499 1
90.5467141139985 43.3906018065003 1
34.5245138532001 60.3963424583717 0
50.2864961189907 49.8045388132306 0
49.5866772163203 59.8089509945327 0
97.6456339600777 68.861572724206 1
32.5772001680931 95.5985476138788 0
74.248691367216 69.8245712265719 1
71.7964620586338 78.4535622451505 1
75.3956114656803 85.7599366733162 1
35.2861128152619 47.0205139472342 0
56.2538174971162 39.2614725105802 0
30.058822446698 49.5929738672369 0
44.6682617248089 66.4500861455891 0
66.5608944724295 41.0920980793697 0
40.4575509837516 97.5351854890994 1
49.0725632190884 51.8832118207397 0
80.27957401467 92.1160608134408 1
66.7467185694404 60.9913940274099 1
32.7228330406032 43.3071730643006 0
64.0393204150601 78.0316880201823 1
72.3464942257992 96.227592967614 1
60.4578857391896 73.0949980975804 1
58.840956217268 75.8584483127904 1
99.8278577969213 72.3692519338389 1
47.2642691084817 88.4758649955978 1
50.4581598028599 75.8098595298246 1
60.4555562927153 42.5084094357222 0
82.2266615778557 42.7198785371646 0
88.9138964166533 69.8037888983547 1
94.834506724302 45.6943068025075 1
67.3192574691753 66.5893531774792 1
57.2387063156986 59.5142819801296 1
80.3667560017127 90.9601478974695 1
68.4685217859111 85.5943071045201 1
42.0754545384731 78.8447860014804 0
75.4777020053391 90.4245389975396 1
78.6354243489802 96.6474271688564 1
52.3480039879411 60.7695052560259 0
94.0943311251679 77.1591050907389 1
90.4485509709636 87.508791764847 1
55.4821611406959 35.5707034722887 0
74.4926924184304 84.8451368493014 1
89.8458067072098 45.3582836109166 1
83.4891627449824 48.3802857972818 1
42.2617008099817 87.1038509402546 1
99.3150088051039 68.7754094720662 1
55.340017560037 64.9319380069486 1
74.7758930009277 89.5298128951328 1];
X = data(:,1:2); % 大写表示矩阵
y = data(:,3); % 小写表示向量
% 初始化参数
[m,n] = size(X);
X = [ones(m,1) X]; % 添加截距项
initial_theta = zeros(n + 1, 1); % 初始化权重参数
iter = 0;
theta = initial_theta;
alpha = 0.001; % 学习率因子
% 梯度下降法优化参数
J_next = 0;
while true
iter = iter + 1;
Z = X * theta;
h = sigmoid(Z); % 调用sigmoid函数计算概率
J = 1 / m * ( - y' * log(h) - (1 - y') * log(1 - h)); % 计算损失函数
grad = 1 / m * (X' * (h - y)); % 梯度
theta = theta - alpha * grad;
if(iter > 1e+10) || abs(J - J_next) < 1e-9
break;
end
J_next = J;
% draw(X(:,2:3),y,theta) % 这是绘图函数,如果不需要可以注释
% pause(0.1)
% hold off
end
disp('三个参数为:')
disp(theta)
disp(['损失函数的值为:',num2str(J)])
draw(X(:,2:3),y,theta)
toc
Sigmoid函数:
function g = sigmoid(Z)
% sigmoid 函数计算
g = zeros(size(Z));
g = 1./(1+exp(-Z));
end
绘图函数:
function draw(X,y,theta)
% 绘制中间过程图
plotData(X,y); % 调用绘图函数
xlabel('score1')
ylabel('score2')
title('梯度下降法分类散点图')
legend('通过','未通过')
% 根据两点确定一条直线,所以只需要找到两个点就可以了
plot_x = [min(X(:,2))-2, max(X(:,2))+2];
plot_y = (-1./theta(3)).*(theta(2).*plot_x + theta(1));
plot(plot_x, plot_y)
axis([30, 100, 30, 100])
end
上面的代码基本上都是基础,中间计算梯度那里用的是矩阵运算。计算结果如下:
2.牛顿法实现逻辑回归
数据仍然是梯度下降法的数据,这里只是把迭代过程改为牛顿法。
主函数:
clc,clear % 清除命令
% 导入数据
tic
data = [];% 梯度下降法的数据
X = data(:,1:2); % 大写表示矩阵
y = data(:,3); % 小写表示向量
% 初始化参数
[m,n] = size(X);
X = [ones(m,1) X]; % 添加截距项
initial_theta = zeros(n + 1, 1); % 初始化权重参数
iter = 0;
theta = initial_theta;
J_next = 0;
%% 牛顿法优化参数
while true
iter = iter + 1;
Z = X * theta;
h = sigmoid(Z); % 获得概率的值
J = 1/m * (-y'*log(h)-(1-y')*log(1-h)); % 这是原函数
grad_1 = 1/m * (X'*(h-y)); % 一阶偏导
grad_2 = h'*(1-h)* X'* X; % 二阶偏导
dx = - grad_2 \ grad_1; % 计算梯度
theta = theta + dx; % 迭代
if (iter>1e+10)||abs(J-J_next) < 1e-9 % 判断条件
break;
end
J_next = J;
end
disp('三个参数为:')
disp(theta)
disp(['损失函数的值为:',num2str(J)])
disp(['迭代次数为:',num2str(iter)])
draw(X(:,2:3),y,theta) % 绘图
toc
绘图函数:
function draw(X,y,theta)
% 绘制中间过程图
plotData(X,y); % 调用绘图函数
xlabel('score1')
ylabel('score2')
title('牛顿法法分类散点图')
legend('通过','未通过')
% 根据两点确定一条直线,所以只需要找到两个点就可以了
plot_x = [min(X(:,2))-2, max(X(:,2))+2];
plot_y = (-1./theta(3)).*(theta(2).*plot_x + theta(1));
plot(plot_x, plot_y)
axis([30, 100, 30, 100])
end
Sigmoid函数:
function g = sigmoid(Z)
% sigmoid 函数计算
g = zeros(size(Z));
g = 1./(1+exp(-Z));
end
计算结果:
牛顿法迭代只用了5.16秒,而梯度下降法迭代使用了近32秒,牛顿法快了近6倍多。迭代次数也要比梯度下降法的迭代次数要少很多,连结果都比梯度下降法要更好。这说明牛顿法相对于梯度下降法来说,收敛效果更好。
总结
逻辑回归是机器学习的一个入门算法,其数学思想较为简单,主要涉及到一个sigmoid函数和似然函数。将极大似然函数转换为极小化损失函数,可采用梯度下降法和牛顿法进行求解,结果表面牛顿法在求解小规模的逻辑回归问题上更有优势。但当样本量增大后,牛顿法的计算量将会很大很大,到时候可能就会用到小批量随机梯度下降法等优化算法了。