【机器学习】逻辑回归实现与应用
介绍
逻辑回归(Logistic Regression),又叫逻辑斯蒂回归,是机器学习中一种十分基础的分类方法,由于算法简单而高效,在实际场景中得到了广泛的应用。本次实验中,我们将探索逻辑回归的原理及算法实现,并使用 scikit-learn 构建逻辑回归分类预测模型。
Sigmoid分布函数
因为今天要解决的是一个二分类问题,所以先介绍此函数
f
(
z
)
=
1
1
+
e
−
z
f(z)=\frac{1}{1+e^{-z}}
f(z)=1+e−z1
python实现:
def sigmoid(z):
sigmoid = 1 / (1 + np.exp(-z))
return sigmoid
这个图像呈现出完美的 S 型(Sigmoid 的含义)。它的取值仅介于
0
0
0 和
1
1
1 之间,且关于
z
=
0
z=0
z=0 轴中心对称。同时当
z
z
z 越大时,
y
y
y 越接近于
1
1
1,而
z
z
z 越小时,
y
y
y 越接近于
0
0
0。如果我们以
0.5
0.5
0.5 为分界点,将
>
0.5
>0.5
>0.5 或
<
0.5
<0.5
<0.5 的值分为两类,这不就是解决
0
−
1
0-1
0−1 二分类问题的完美选择嘛。
逻辑回归模型
如果一组连续随机变量符合 Sigmoid 函数样本分布,就称作为逻辑分布。逻辑分布是概率论中的定理,是一种连续型的概率分布。
在逻辑回归中,定义:
z
i
=
w
0
x
0
+
w
1
x
1
+
⋯
+
w
i
x
i
=
w
T
x
z_{i} = {w_0}{x_0} + {w_1}{x_1} + \cdots + {w_i}{x_i} = {w^T}x
zi=w0x0+w1x1+⋯+wixi=wTx
f
(
z
i
)
=
1
1
+
e
−
z
i
f(z_{i})=\frac{1}{1+e^{-z_{i}}}
f(zi)=1+e−zi1
即:
h
w
(
x
)
=
f
(
w
T
x
)
=
1
1
+
e
−
w
T
x
h_{w}(x) = f({w^T}x)=\frac{1}{1+e^{-w^Tx}}
hw(x)=f(wTx)=1+e−wTx1
由于目标值
y
y
y 只有 0 和 1 两个值,那么如果记
y
=
1
y=1
y=1 的概率为
h
w
(
x
)
h_{w}(x)
hw(x),则此时
y
=
0
y=0
y=0 的概率为
1
−
h
w
(
x
)
1-h_{w}(x)
1−hw(x)。那么,我们可以记作逻辑回归模型条件概率分布:
P
(
Y
=
y
∣
x
)
=
{
h
w
(
x
)
,
y
=
1
1
−
h
w
(
x
)
,
y
=
0
P(Y=y | x)=\left\{\begin{array}{rlrl}{h_{w}(x)} & {, y=1} \\ {1-h_{w}(x)} & {, y=0}\end{array}\right.
P(Y=y∣x)={hw(x)1−hw(x),y=1,y=0
其可等价写为似然函数:
P
(
y
∣
x
;
w
)
=
(
h
w
(
x
)
)
y
(
1
−
h
w
(
x
)
)
1
−
y
P(y|x ; w)=\left(h_{w}(x)\right)^{y}\left(1-h_{w}(x)\right)^{1-y}
P(y∣x;w)=(hw(x))y(1−hw(x))1−y
对于
i
i
i 个样本的总概率而言实际上可以看作单样本概率的乘积,记为
L
(
w
)
L(w)
L(w):
L
(
w
)
=
∏
i
=
1
m
(
h
w
(
x
(
i
)
)
)
y
(
i
)
(
1
−
h
w
(
x
(
i
)
)
)
1
−
y
(
i
)
L(w) =\prod_{i=1}^{m}\left(h_{w}\left(x^{(i)}\right)\right)^{y^{(i)}}\left(1-h_{w}\left(x^{(i)}\right)\right)^{1-y^{(i)}}
L(w)=i=1∏m(hw(x(i)))y(i)(1−hw(x(i)))1−y(i)
由于连乘表示起来非常复杂,我们应用数学技巧,即两边取对数将连乘转换为连加的形式,即:
log
L
(
w
)
=
∑
i
=
1
m
[
y
(
i
)
log
h
w
(
x
(
i
)
)
+
(
1
−
y
(
i
)
)
log
(
1
−
h
w
(
x
(
i
)
)
)
]
(1)
\log L(w)=\sum_{i=1}^{m} \left [ y^{(i)} \log h_{w}\left(x^{(i)}\right)+\left(1-y^{(i)}\right) \log \left(1-h_{w}\left(x^{(i)}\right)\right)\right ] \tag{1}
logL(w)=i=1∑m[y(i)loghw(x(i))+(1−y(i))log(1−hw(x(i)))](1)
对数损失函数
函数(1)衡量了事件发生的总概率。根据最大似然估计原理,只需要通过对 𝐿(𝑤) 求最大值,即得到 𝑤 的估计值。而在机器学习问题中,我们需要一个损失函数,并通过求其最小值来进行参数优化。所以,对数似然函数取负数就可以被作为逻辑回归的对数损失函数:
J
(
w
)
=
−
1
m
∑
i
=
1
m
[
y
(
i
)
log
h
w
(
x
(
i
)
)
+
(
1
−
y
(
i
)
)
log
(
1
−
h
w
(
x
(
i
)
)
)
]
J(w) =- \frac{1}{m} \sum_{i=1}^{m} \left [ y^{(i)} \log h_{w}\left(x^{(i)}\right)+\left(1-y^{(i)}\right) \log \left(1-h_{w}\left(x^{(i)}\right)\right)\right ]
J(w)=−m1i=1∑m[y(i)loghw(x(i))+(1−y(i))log(1−hw(x(i)))]
python实现:
def loss(h, y):
loss = (-y * np.log(h) - (1 - y) * np.log(1 - h)).mean()
return loss
梯度下降法
之前我已经写过一期梯度下降法相关的blog,名为【机器学习】梯度下降法求解线性回归参数,但是其中对于数学原理的介绍并不是很详尽,这次我将结合上述本次公式对梯度下降法再次讲解,其实梯度的概念是高等数学中的内容,很明显重点就是求导,因为我们想知道的是如何下降最快,这样只需要沿着上升最快的反方向走就可以了,也就是沿着梯度的反向走,按照步长一点一点走,不然就有可能跳过最小值。
梯度是一个向量,它表示某一函数在该点处的方向导数沿着该方向取得最大值,即函数在该点处沿着该方向(此梯度的方向)变化最快,变化率最大(为该梯度的模)。简而言之,对于一元函数而言,梯度就是指在某一点的导数。而对于多元函数而言,梯度就是指在某一点的偏导数组成的向量。
我们先针对公式(1)化简:
J
(
w
)
=
−
1
m
∑
i
=
1
m
[
y
(
i
)
(
−
log
(
1
+
e
−
w
T
x
(
i
)
)
)
+
(
1
−
y
(
i
)
)
(
−
w
T
x
(
i
)
−
log
(
1
+
e
−
w
T
x
(
i
)
)
)
]
=
−
1
m
∑
i
=
1
m
[
−
y
(
i
)
log
(
1
+
e
−
w
T
x
(
i
)
)
−
w
T
x
(
i
)
−
log
(
1
+
e
−
w
T
x
(
i
)
)
+
w
T
x
(
i
)
y
(
i
)
+
y
(
i
)
log
(
1
+
e
−
w
T
x
(
i
)
)
]
=
−
1
m
∑
i
=
1
m
[
−
w
T
x
(
i
)
−
log
(
1
+
e
−
w
T
x
(
i
)
)
+
w
T
x
(
i
)
y
(
i
)
]
=
−
1
m
∑
i
=
1
m
[
−
log
(
e
w
T
x
(
i
)
)
−
log
(
1
+
e
−
w
T
x
(
i
)
)
+
w
T
x
(
i
)
y
(
i
)
]
=
−
1
m
∑
i
=
1
m
[
−
log
(
e
w
T
x
(
i
)
(
1
+
e
−
w
T
x
(
i
)
)
)
+
w
T
x
(
i
)
y
(
i
)
]
=
−
1
m
∑
i
=
1
m
[
−
log
(
e
w
T
x
(
i
)
+
1
)
+
w
T
x
(
i
)
y
(
i
)
]
(2)
\begin{aligned} J(w) & = -\frac{1}{m} \sum_{i=1}^m \left [ y^{(i)} (-\log(1+e^{-w^Tx^{(i)}})) + (1-y^{(i)})(-w^Tx^{(i)}-\log(1+e^{-w^Tx^{(i)}}))\right ] \\ & = -\frac{1}{m} \sum_{i=1}^m \left [ - y^{(i)}\log(1+e^{-w^Tx^{(i)}}) - w^Tx^{(i)}-\log(1+e^{-w^Tx^{(i)}})+w^Tx^{(i)}y^{(i)}+y^{(i)}\log(1+e^{-w^Tx^{(i)}})\right ]\\ & = -\frac{1}{m} \sum_{i=1}^m \left [- w^Tx^{(i)}-\log(1+e^{-w^Tx^{(i)}})+w^Tx^{(i)}y^{(i)} \right ]\\ & = -\frac{1}{m} \sum_{i=1}^m \left [-\log(e^{w^Tx^{(i)}})-\log(1+e^{-w^Tx^{(i)}})+w^Tx^{(i)}y^{(i)} \right ]\\ & = -\frac{1}{m} \sum_{i=1}^m \left [-\log(e^{w^Tx^{(i)}}(1+e^{-w^Tx^{(i)}}))+w^Tx^{(i)}y^{(i)} \right ]\\ & = -\frac{1}{m} \sum_{i=1}^m \left [-\log(e^{w^Tx^{(i)}}+1)+w^Tx^{(i)}y^{(i)} \right ]\tag{2}\\ \end{aligned}
J(w)=−m1i=1∑m[y(i)(−log(1+e−wTx(i)))+(1−y(i))(−wTx(i)−log(1+e−wTx(i)))]=−m1i=1∑m[−y(i)log(1+e−wTx(i))−wTx(i)−log(1+e−wTx(i))+wTx(i)y(i)+y(i)log(1+e−wTx(i))]=−m1i=1∑m[−wTx(i)−log(1+e−wTx(i))+wTx(i)y(i)]=−m1i=1∑m[−log(ewTx(i))−log(1+e−wTx(i))+wTx(i)y(i)]=−m1i=1∑m[−log(ewTx(i)(1+e−wTx(i)))+wTx(i)y(i)]=−m1i=1∑m[−log(ewTx(i)+1)+wTx(i)y(i)](2)
接下来,我们针对公式(2)求导:
∂
J
∂
w
=
−
1
m
∑
i
=
1
m
[
−
x
(
i
)
e
w
T
x
(
i
)
e
w
T
x
(
i
)
+
1
+
x
(
i
)
y
(
i
)
]
=
−
1
m
∑
i
=
1
m
[
−
e
w
T
x
(
i
)
e
w
T
x
(
i
)
+
1
+
y
(
i
)
]
x
(
i
)
=
−
1
m
∑
i
=
1
m
[
−
e
−
w
T
x
(
i
)
e
w
T
x
(
i
)
e
−
w
T
x
(
i
)
(
e
w
T
x
(
i
)
+
1
)
+
y
(
i
)
]
x
(
i
)
=
−
1
m
∑
i
=
1
m
[
−
1
1
+
e
−
w
T
x
(
i
)
+
y
(
i
)
]
x
(
i
)
\begin{aligned} \frac{\partial{J}}{\partial{w}} & = -\frac{1}{m} \sum_{i=1}^m \left [ - \frac{x^{(i)}e^{w^Tx^{(i)}}}{e^{w^Tx^{(i)}}+1}+x^{(i)}y^{(i)}\right ]\\ & = -\frac{1}{m} \sum_{i=1}^m \left [ - \frac{e^{w^Tx^{(i)}}}{e^{w^Tx^{(i)}}+1}+y^{(i)}\right ]x^{(i)}\\ & = -\frac{1}{m} \sum_{i=1}^m \left [ - \frac{e^{-w^Tx^{(i)}}e^{w^Tx^{(i)}}}{e^{-w^Tx^{(i)}}(e^{w^Tx^{(i)}}+1)}+y^{(i)}\right ]x^{(i)}\\ & = -\frac{1}{m} \sum_{i=1}^m \left [ - \frac{1}{1+e^{-w^Tx^{(i)}}}+y^{(i)}\right ]x^{(i)}\\ \end{aligned}
∂w∂J=−m1i=1∑m[−ewTx(i)+1x(i)ewTx(i)+x(i)y(i)]=−m1i=1∑m[−ewTx(i)+1ewTx(i)+y(i)]x(i)=−m1i=1∑m[−e−wTx(i)(ewTx(i)+1)e−wTx(i)ewTx(i)+y(i)]x(i)=−m1i=1∑m[−1+e−wTx(i)1+y(i)]x(i)
即:
∂
J
∂
w
=
−
1
m
∑
i
=
1
m
[
−
h
w
(
x
(
i
)
)
+
y
(
i
)
]
x
(
i
)
\frac{\partial{J}}{\partial{w}}= -\frac{1}{m} \sum_{i=1}^m \left [ - h_w(x^{(i)})+y^{(i)}\right ]x^{(i)}
∂w∂J=−m1i=1∑m[−hw(x(i))+y(i)]x(i)
向量的形式表达为:
∂
J
∂
w
=
1
m
x
T
(
h
w
(
x
)
−
y
)
\frac{\partial{J}}{\partial{w}} = \frac{1}{m}x^T(h_{w}(x)-y)
∂w∂J=m1xT(hw(x)−y)
当我们得到梯度的方向,然后乘以一个常数
α
\alpha
α,就可以得到每次梯度下降的步长(上图箭头的长度)。最后,通过多次迭代,找到梯度变化很小的点,也就对应着损失函数的极小值了。其中,常数
α
\alpha
α 往往也被称之为学习率 Learning Rate。执行权重更新的过程为:
w
←
w
−
α
∂
J
∂
w
w \leftarrow w - \alpha \frac{\partial{J}}{\partial{w}}
w←w−α∂w∂J
python实现:
def gradient(X, h, y):
gradient = np.dot(X.T, (h - y)) / y.shape[0]
return gradient
逻辑回归实现
加载数据
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
df = pd.read_csv(
"https://labfile.oss.aliyuncs.com/courses/1081/course-8-data.csv", header=0) # 加载数据集
df.head() # 预览前 5 行数据
函数代码汇总
def sigmoid(z):
# Sigmoid 分布函数
sigmoid = 1 / (1 + np.exp(-z))
return sigmoid
def loss(h, y):
# 损失函数
loss = (-y * np.log(h) - (1 - y) * np.log(1 - h)).mean()
return loss
def gradient(X, h, y):
# 梯度计算
gradient = np.dot(X.T, (h - y)) / y.shape[0]
return gradient
def Logistic_Regression(x, y, lr, num_iter):
# 逻辑回归过程
intercept = np.ones((x.shape[0], 1)) # 初始化截距为 1
x = np.concatenate((intercept, x), axis=1)
w = np.zeros(x.shape[1]) # 初始化参数为 0
for i in range(num_iter): # 梯度下降迭代
z = np.dot(x, w) # 线性函数
h = sigmoid(z) # sigmoid 函数
g = gradient(x, h, y) # 计算梯度
w -= lr * g # 通过学习率 lr 计算步长并执行梯度下降
l = loss(h, y) # 计算损失函数值
return l, w # 返回迭代后的梯度和参数
逻辑回归
x = df[['X0', 'X1']].values
y = df['Y'].values
lr = 0.01 # 学习率
num_iter = 30000 # 迭代次数
# 训练
L = Logistic_Regression(x, y, lr, num_iter)
L
逻辑回归 scikit-learn 实现
LogisticRegression(penalty='l2', dual=False, tol=0.0001, C=1.0, fit_intercept=True, intercept_scaling=1, class_weight=None, random_state=None, solver='liblinear', max_iter=100, multi_class='ovr', verbose=0, warm_start=False, n_jobs=1)介绍其中几个常用的参数,其余使用默认即可:
penalty
: 惩罚项,默认为 L 2 L_{2} L2 范数。dual
: 对偶化,默认为 False。tol
: 数据解算精度。fit_intercept
: 默认为 True,计算截距项。random_state
: 随机数发生器。max_iter
: 最大迭代次数,默认为 100。另外,
solver
参数用于指定求解损失函数的方法。默认为liblinear
(0.22 开始默认为
lbfgs
),适合于小数据集。除此之外,还有适合多分类问题的newton-cg
,sag
,saga
和lbfgs
求解器。这些方法都来自于一些学术论文,有兴趣可以自行搜索了解。
from sklearn.linear_model import LogisticRegression
model = LogisticRegression(tol=0.001, max_iter=10000, solver='liblinear') # 设置数据解算精度和迭代次数
model.fit(x, y)
model.coef_, model.intercept_