微信公众号 slandarer随笔 中
这篇文章正在抽奖赠书
欢迎过去参与抽奖，文末有该书更详细介绍：
请在微信打开抽奖地址
抽奖推送链接：抽奖地址
之前的最小二乘法的两种解读那篇文章，我们拟合了多项式，拟合了线性多元函数，我们的函数映射结果是数值，但我们想要的得到的结果如果不是数值，而是(是/否)(TRUE/FALSE)应该怎么做?此即逻辑回归。
代价函数及其梯度(Cost function and gradient)
首先来看要拟合的函数，下面是我们之前拟合的函数的形式，很明显拟合结果为一个数值。
h
θ
(
x
)
=
θ
T
x
=
∑
i
=
0
n
θ
i
x
i
=
θ
0
+
θ
1
x
1
+
?
+
θ
n
x
n
h_{\theta}(x)=\theta^Tx=\sum_{i=0}^{n} \theta_{i} x_{i}=\theta_{0}+\theta_{1} x_{1}+\cdots+\theta_{n} x_{n}
hθ?(x)=θTx=i=0∑n?θi?xi?=θ0?+θ1?x1?+?+θn?xn?
我们想要将一个范围
(
?
∞
,
∞
)
(-\infty,\infty)
(?∞,∞)的数值映射到(是/否)(1/0)我们非常自然的能够想到sigmoid函数，这是一个能将数值从
(
?
∞
,
∞
)
(-\infty,\infty)
(?∞,∞)映射到
(
0
,
1
)
(0,1)
(0,1)的函数：
g
(
z
)
=
1
1
+
e
?
z
g(z)=\frac{1}{1+e^{-z}}
g(z)=1+e?z1?
我们直接把sigmoid函数往原本的拟合函数外面一套，拟合的函数不就能把变量很顺滑的映射到
(
0
,
1
)
(0,1)
(0,1)了嘛，然后我们认为大于0.5就代表是，小于0.5就代表否，美滋滋：
h
θ
(
x
)
=
g
(
θ
T
x
)
h_{\theta}(x)=g(\theta^Tx)
hθ?(x)=g(θTx)
我们再来看代价函数的改变，代价函数我们并没用平常的平方损失函数，而是用了如下的一个形式：
J
(
θ
)
=
1
m
∑
i
=
1
m
[
?
y
(
i
)
log
?
(
h
θ
(
x
(
i
)
)
)
?
(
1
?
y
(
i
)
)
log
?
(
1
?
h
θ
(
x
(
i
)
)
)
]
J(\theta)=\frac{1}{m} \sum_{i=1}^{m}\left[-y^{(i)} \log \left(h_{\theta}\left(x^{(i)}\right)\right)-\left(1-y^{(i)}\right) \log \left(1-h_{\theta}\left(x^{(i)}\right)\right)\right]
J(θ)=m1?i=1∑m?[?y(i)log(hθ?(x(i)))?(1?y(i))log(1?hθ?(x(i)))]
别看这个形式有些复杂，但其实并没有那么可怕，首先它将m个数据组都带入了后面的式子求了一个均值。而对于单独一组数据，我们的代价函数的形式是下面这样的：
J
(
θ
)
=
?
y
log
?
(
h
θ
(
x
)
)
?
(
1
?
y
)
log
?
(
1
?
h
θ
(
x
)
)
J(\theta)=-y\log \left(h_{\theta}\left(x\right)\right)-\left(1-y\right) \log \left(1-h_{\theta}\left(x\right)\right)
J(θ)=?ylog(hθ?(x))?(1?y)log(1?hθ?(x))
我们的y值取值只有(1/0)，x很难真正取到正负无穷，也就是说
h
θ
(
x
)
h_{\theta}(x)
hθ?(x)范围为(0,1)并不会真正取到两端点，这时候代价函数可以看作一个条件函数：
J
(
θ
)
=
{
?
log
?
(
h
θ
(
x
)
)
y
=
1
?
log
?
(
1
?
h
θ
(
x
)
)
y
=
0
J(\theta)=\begin{cases} -\log \left(h_{\theta}\left(x\right)\right) &y=1\\ -\log \left(1-h_{\theta}\left(x\right)\right) &y=0 \end{cases}
J(θ)={?log(hθ?(x))?log(1?hθ?(x))?y=1y=0?
y=1时，我们考察的是
h
θ
(
x
)
h_{\theta}(x)
hθ?(x)与1的接近程度，越接近1，
?
log
?
(
h
θ
(
x
)
)
-\log(h_{\theta}(x))
?log(hθ?(x))就越趋紧于0，y=0时，我们考察的是
1
?
h
θ
(
x
)
1-h_{\theta}(x)
1?hθ?(x)与1的接近程度，越接近1，
?
log
?
(
1
?
h
θ
(
x
)
)
-\log(1-h_{\theta}(x))
?log(1?hθ?(x))就越趋紧于0，总而言之，这是一个越接近真实值，计算值越接近0，越远离真实值，计算值越接近正无穷的非常巧妙的函数。
我们要求代价函数的极小值点，当然还是要让偏导等于0，这里我们对 J ( θ ) J(\theta) J(θ)求一下偏导，这玩意没那么好求，因此我们就一层一层来呗，先研究对sigmoid函数求导，再研究对 h θ ( x ) h_{\theta}(x) hθ?(x)求偏导然后链式求导法则走起来：
第三层求导：
g
(
z
)
=
1
1
+
e
?
z
g
′
(
z
)
=
e
?
z
(
1
+
e
?
z
)
2
=
1
+
e
?
z
?
1
(
1
+
e
?
z
)
2
=
1
1
+
e
?
z
(
1
+
e
?
z
?
1
1
+
e
?
z
)
=
g
(
z
)
(
1
?
g
(
z
)
)
\begin{aligned} g(z)&=\frac{1}{1+e^{-z}}\\ g'(z)&=\frac{e^{-z}}{(1+e^{-z})^2}=\frac{1+e^{-z}-1}{(1+e^{-z})^2}\\ &=\frac{1}{1+e^{-z}}(\frac{1+e^{-z}-1}{1+e^{-z}})\\ &=g(z)(1-g(z)) \end{aligned}
g(z)g′(z)?=1+e?z1?=(1+e?z)2e?z?=(1+e?z)21+e?z?1?=1+e?z1?(1+e?z1+e?z?1?)=g(z)(1?g(z))?
第二层求导：
h
θ
(
x
)
=
g
(
θ
T
x
)
=
1
1
+
e
?
θ
T
x
=
1
1
+
e
?
(
θ
0
+
θ
1
x
1
+
?
+
θ
n
x
n
)
?
h
θ
(
x
)
?
θ
j
=
1
1
+
e
?
(
θ
0
+
θ
1
x
1
+
?
+
θ
n
x
n
)
(
1
?
1
1
+
e
?
(
θ
0
+
θ
1
x
1
+
?
+
θ
n
x
n
)
)
x
j
=
h
θ
(
x
)
(
1
?
h
θ
(
x
)
)
x
j
?
l
o
g
(
h
θ
(
x
)
)
?
θ
j
=
h
θ
(
x
)
(
1
?
h
θ
(
x
)
)
x
j
h
θ
(
x
)
=
(
1
?
h
θ
(
x
)
)
x
j
?
l
o
g
(
1
?
h
θ
(
x
)
)
?
θ
j
=
?
h
θ
(
x
)
(
1
?
h
θ
(
x
)
)
x
j
1
?
h
θ
(
x
)
=
?
h
θ
(
x
)
x
j
\begin{aligned} h_{\theta}(x)&=g(\theta^Tx)=\frac{1}{1+e^{-\theta^Tx}}=\frac{1}{1+e^{-(\theta_{0}+\theta_{1} x_{1}+\cdots+\theta_{n} x_{n})}}\\ \\ \frac{\partial h_{\theta}(x)}{\partial\theta_{j}}&= \frac{1}{1+e^{-(\theta_{0}+\theta_{1} x_{1}+\cdots+\theta_{n} x_{n})}}\left(1-\frac{1}{1+e^{-(\theta_{0}+\theta_{1} x_{1}+\cdots+\theta_{n} x_{n})}}\right)x_j\\ &=h_{\theta}(x)(1-h_{\theta}(x))x_j\\ \\ \frac{\partial log(h_{\theta}(x))}{\partial\theta_{j}}&= \frac{h_{\theta}(x)(1-h_{\theta}(x))x_j}{h_{\theta}(x)}\\ &=(1-h_{\theta}(x))x_j\\ \\ \frac{\partial log(1-h_{\theta}(x))}{\partial\theta_{j}}&= \frac{-h_{\theta}(x)(1-h_{\theta}(x))x_j}{1-h_{\theta}(x)}\\ &=-h_{\theta}(x)x_j \end{aligned}
hθ?(x)?θj??hθ?(x)??θj??log(hθ?(x))??θj??log(1?hθ?(x))??=g(θTx)=1+e?θTx1?=1+e?(θ0?+θ1?x1?+?+θn?xn?)1?=1+e?(θ0?+θ1?x1?+?+θn?xn?)1?(1?1+e?(θ0?+θ1?x1?+?+θn?xn?)1?)xj?=hθ?(x)(1?hθ?(x))xj?=hθ?(x)hθ?(x)(1?hθ?(x))xj??=(1?hθ?(x))xj?=1?hθ?(x)?hθ?(x)(1?hθ?(x))xj??=?hθ?(x)xj??
第一层求导：
先看看对于一组数据：
J
(
θ
)
=
?
y
log
?
(
h
θ
(
x
)
)
?
(
1
?
y
)
log
?
(
1
?
h
θ
(
x
)
)
?
J
(
θ
)
?
θ
j
=
?
y
(
1
?
h
θ
(
x
)
)
x
j
+
(
1
?
y
)
h
θ
(
x
)
x
j
=
?
y
x
j
+
y
h
θ
(
x
)
x
j
+
h
θ
(
x
)
x
j
?
y
h
θ
(
x
)
x
j
=
(
h
θ
(
x
)
?
y
)
x
j
\begin{aligned} J(\theta)&=-y\log(h_{\theta}(x))-(1-y)\log(1-h_{\theta}(x))\\ \frac{\partial J(\theta)}{\partial\theta_{j}}&= -y(1-h_{\theta}(x))x_j+(1-y)h_{\theta}(x)x_j\\ &=-yx_j+yh_{\theta}(x)x_j+h_{\theta}(x)x_j-yh_{\theta}(x)x_j\\ &=(h_{\theta}(x)-y)x_j \end{aligned}
J(θ)?θj??J(θ)??=?ylog(hθ?(x))?(1?y)log(1?hθ?(x))=?y(1?hθ?(x))xj?+(1?y)hθ?(x)xj?=?yxj?+yhθ?(x)xj?+hθ?(x)xj??yhθ?(x)xj?=(hθ?(x)?y)xj??
对于m组数据则有：
?
J
(
θ
)
?
θ
j
=
1
m
∑
i
=
1
m
(
h
θ
(
x
(
i
)
)
?
y
(
i
)
)
x
j
(
i
)
\frac{\partial J(\theta)}{\partial \theta_{j}}=\frac{1}{m} \sum_{i=1}^{m}\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right) x_{j}^{(i)}
?θj??J(θ)?=m1?i=1∑m?(hθ?(x(i))?y(i))xj(i)?
特征映射(Feature mapping)
其次再看特征映射，对于如下情况我们假设横坐标为
x
1
x_1
x1?，纵坐标为
x
2
x_2
x2?:
很明显分类结果不光和 x 1 , x 2 x_1,x_2 x1?,x2?有关，还和交叉项有关，其次我们想到之前最小二乘法那篇时，我们使用了 1 , x , x 2 , x 3 , … 1,x,x^2,x^3,\dots 1,x,x2,x3,…一系列项来拟合函数，在这里我们也很自然的能想到用 1 , x 1 , x 2 , x 1 2 , x 1 x 2 , x 2 2 , x 1 3 , x 1 2 x 2 , … 1,x_1,x_2,x_1^2,x_1x_2,x_2^2,x_1^3,x_1^2x_2,\dots 1,x1?,x2?,x12?,x1?x2?,x22?,x13?,x12?x2?,…等一系列项来拟合函数，这里我们用了的项的最高次为6次，因而一共有28项：
m a p F e a t u r e ( x ) = ( ∣ ∣ ∣ ∣ ∣ ∣ ∣ … 1 x 1 x 2 x 1 2 x 1 x 2 x 2 2 x 1 3 … ∣ ∣ ∣ ∣ ∣ ∣ ∣ … ) mapFeature(x)=\begin{pmatrix} \mid & \mid & \mid & \mid & \mid & \mid & \mid &\dots \\ 1& x_1& x_2& x_1^2& x_1x_2& x_2^2& x_1^3&\dots \\ \mid & \mid & \mid & \mid & \mid & \mid & \mid &\dots \end{pmatrix} mapFeature(x)=???∣1∣?∣x1?∣?∣x2?∣?∣x12?∣?∣x1?x2?∣?∣x22?∣?∣x13?∣?………????
正则化逻辑回归(Regularized logistic regression)
讲明白啥是逻辑回归，这里再说明一下正则化的作用，我们要对逻辑回归进行正则化首先是要防止其过拟合和欠拟合(主要是过拟合)，就是防止图中左图及右图的情况：
什么时候容易出现左图情况？就是拟合所用的参数数量比较少，或者说大部分系数数值比较小的情况，而出现右图这种情况一般是参数数量比较多，或者说部分没那么重要的系数数值比较大的情况，因此我们想要拟合结果合理，我们有两种途径，一是调整参数数量，另一个便是将系数数值限定在一定范围内，这里正则化就是用到的第二种方法，我们将代价函数改写为：
J
(
θ
)
=
1
m
∑
i
=
1
m
[
?
y
(
i
)
log
?
(
h
θ
(
x
(
i
)
)
)
?
(
1
?
y
(
i
)
)
log
?
(
1
?
h
θ
(
x
(
i
)
)
)
]
+
λ
2
m
∑
j
=
1
n
θ
j
2
J(\theta)=\frac{1}{m} \sum_{i=1}^{m}\left[-y^{(i)} \log \left(h_{\theta}\left(x^{(i)}\right)\right)-\left(1-y^{(i)}\right) \log \left(1-h_{\theta}\left(x^{(i)}\right)\right)\right]+\frac{\lambda}{2 m} \sum_{j=1}^{n} \theta_{j}^{2}
J(θ)=m1?i=1∑m?[?y(i)log(hθ?(x(i)))?(1?y(i))log(1?hθ?(x(i)))]+2mλ?j=1∑n?θj2?
这样就能起到限制
θ
\theta
θ大小的作用，需要注意的是
λ
\lambda
λ过大容易让
θ
\theta
θ过小从而出现欠拟合，
λ
\lambda
λ过小甚至为0容易让
θ
\theta
θ过大从而出现过拟合，因此调节合适的
λ
\lambda
λ大小非常重要。相应的偏导也变为：
?
J
(
θ
)
?
θ
j
=
(
1
m
∑
i
=
1
m
(
h
θ
(
x
(
i
)
)
?
y
(
i
)
)
x
j
(
i
)
)
+
λ
m
θ
j
\frac{\partial J(\theta)}{\partial \theta_{j}}=\left(\frac{1}{m} \sum_{i=1}^{m}\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right) x_{j}^{(i)}\right)+\frac{\lambda}{m} \theta_{j}
?θj??J(θ)?=(m1?i=1∑m?(hθ?(x(i))?y(i))xj(i)?)+mλ?θj?
代码部分
数据来源于吴恩达的Machine Learning | Programming Exercise 2: Logistic Regression,这个练习2的完整文件及数据也会和完整代码一块放到文末提供的压缩包内。
数据是某一制造厂微芯片的两次测试结果及合格情况，数据第一第二列为两次测试结果，第三列为合格情况，合格记为1，不合格记为0。
导入数据并绘图
% 数据导入
oriData=readmatrix('ex2data2.txt');
X=oriData(:,[1,2]);
y=oriData(:,3);
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% 绘制图像
ax=gca;
hold(ax,'on')
plot(X(y==0,1),X(y==0,2),'ko','MarkerFaceColor','y','MarkerSize',7)
plot(X(y==1,1),X(y==1,2),'k+','LineWidth',2,'MarkerSize',7)
% 修饰一下
grid(ax,'on')
ax.FontName='cambria';
ax.XColor=[1,1,1].*.3;
ax.YColor=[1,1,1].*.3;
ax.ZColor=[1,1,1].*.3;
ax.LineWidth=1.5;
ax.GridLineStyle='--';
ax.DataAspectRatio=[1,1,1];
数据特征映射
% 将X变量进行扩充为X_exp
% X_exp第一列为数值1
% degree为X1和X2的次数和(阶数)，例如X1*X^2的degree为3
% 我们认为y的取值不光与X1，X2有关，还和其次方项及交叉项有关
% 该过程就是将原本两项[X1,X2],扩充为包含常数项在内的28项
% [X1,X2] => [1,X1,X2,X1^2,X1X2,X2^2,X1^3,X1^2*X2,X1*X2^2,X2^3,......]
degree=6;
X_exp=ones(size(X(:,1)));
for i=1:degree
for j=0:i
X_exp(:,end+1)=(X(:,1).^(i-j)).*(X(:,2).^j);
end
end
构造代价函数及sigmoid函数
% sigmoid函数
function g=sigmoid(z)
g=1./(1+exp(-z));
end
% 梯度及代价值计算函数
function [J,grad]=costFunctionReg(theta,X,y,lambda)
m=length(y);
h_theta=sigmoid(X*theta);
J=sum(-y.*log(h_theta)-(1-y).*log(1-h_theta))/m+sum(theta(2:end).^2)*lambda/2/m;
grad=(sum((h_theta-y).*X,1)./m)'+[0;theta(2:end).*lambda./m];
end
求解
% 设置初始参数和梯度下降率
initial_theta=zeros(size(X_exp,2),1);
lambda=1;
% 设置条件并用自带函数寻找最优参数
options=optimset('GradObj','on','MaxIter',400);
[theta,J,exit_flag]= ...
fminunc(@(t)(costFunctionReg(t,X_exp,y,lambda)),initial_theta,options);
绘制分界曲线
% 绘制分界曲线
degree=6;
X_sym=sym('1');
for i=1:degree
for j=0:i
X_sym(:,end+1)=(sym('x').^(i-j)).*(sym('y').^j);
end
end
Func=matlabFunction(1./(1+exp(X_sym*theta))-.5);
fimplicit(Func,'LineWidth',2,'Color',[169,64,71]./255)
计算正确率
% 将根据梯度下降法得出的预测结果与实际情况对比，并计算正确率
p=sigmoid(X_exp*theta)>.5;
ac=mean(double(p == y))*100;
fprintf('准确率为: %f\n', ac);
完整代码
%@author:slandarer
% 数据导入
oriData=readmatrix('ex2data2.txt');
X=oriData(:,[1,2]);
y=oriData(:,3);
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% 绘制图像
ax=gca;
hold(ax,'on')
plot(X(y==0,1),X(y==0,2),'ko','MarkerFaceColor','y','MarkerSize',7)
plot(X(y==1,1),X(y==1,2),'k+','LineWidth',2,'MarkerSize',7)
% 修饰一下
grid(ax,'on')
ax.FontName='cambria';
ax.XColor=[1,1,1].*.3;
ax.YColor=[1,1,1].*.3;
ax.ZColor=[1,1,1].*.3;
ax.LineWidth=1.5;
ax.GridLineStyle='--';
ax.DataAspectRatio=[1,1,1];
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% 将X变量进行扩充为X_exp
% X_exp第一列为数值1
% degree为X1和X2的次数和(阶数)，例如X1*X^2的degree为3
% 我们认为y的取值不光与X1，X2有关，还和其次方项及交叉项有关
% 该过程就是将原本两项[X1,X2],扩充为包含常数项在内的28项
% [X1,X2] => [1,X1,X2,X1^2,X1X2,X2^2,X1^3,X1^2*X2,X1*X2^2,X2^3,......]
degree=6;
X_exp=ones(size(X(:,1)));
for i=1:degree
for j=0:i
X_exp(:,end+1)=(X(:,1).^(i-j)).*(X(:,2).^j);
end
end
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% 设置初始参数和梯度下降率
initial_theta=zeros(size(X_exp,2),1);
lambda=1;
% 设置条件并用自带函数寻找最优参数
options=optimset('GradObj','on','MaxIter',400);
[theta,J,exit_flag]= ...
fminunc(@(t)(costFunctionReg(t,X_exp,y,lambda)),initial_theta,options);
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% 绘制分界曲线
degree=6;
X_sym=sym('1');
for i=1:degree
for j=0:i
X_sym(:,end+1)=(sym('x').^(i-j)).*(sym('y').^j);
end
end
Func=matlabFunction(1./(1+exp(X_sym*theta))-.5);
fimplicit(Func,'LineWidth',2,'Color',[169,64,71]./255)
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% 将根据梯度下降法得出的预测结果与实际情况对比，并计算正确率
p=sigmoid(X_exp*theta)>.5;
ac=mean(double(p == y))*100;
fprintf('准确率为: %f\n', ac);
% =========================================================================
% sigmoid函数
function g=sigmoid(z)
g=1./(1+exp(-z));
end
% 梯度及代价值计算函数
function [J,grad]=costFunctionReg(theta,X,y,lambda)
m=length(y);
h_theta=sigmoid(X*theta);
J=sum(-y.*log(h_theta)-(1-y).*log(1-h_theta))/m+sum(theta(2:end).^2)*lambda/2/m;
grad=(sum((h_theta-y).*X,1)./m)'+[0;theta(2:end).*lambda./m];
end
完整文件压缩包
链接： https://pan.baidu.com/s/1YIWxRx1qShsbYytjv-yrcA?pwd=slan
提取码： slan
赠书介绍
本书以简单的组合优化问题作为MATLAB智能优化算法实战应用的切入点，逐步深入使用MATLAB编写更复杂的智能优化算法和求解更复杂的组合优化问题，通过9个常见的组合优化问题、5个经典的智能优化算法及4个新颖的智能优化算法，让读者逐渐理解智能优化算法的实际求解过程、算法设计思路及代码编写思路。
目录
（1）遗传算法求解0-1背包问题；
（2）变邻域搜索算法求解旅行商问题
（3）大规模邻域搜索算法求解旅行商问题
（4）灰狼优化算法求解多旅行商问题
（5）蚁群算法求解容量受限的车辆路径问题
（6）模拟退火算法求解同时取送货的车辆路径问题
（7）遗传算法求解带时间窗的车辆路径问题
（8）萤火虫算法求解订单分批问题
（9）头脑风暴优化算法求解带时间窗和同时取送货的车辆路径问题
（10）鲸鱼优化算法求解开放式车辆路径问题。
抽奖地址=〉抽奖地址
微信公众号地址=〉
当然即使没有抽中也可以在如下链接进行购买：https://item.jd.com/13422442.html