在感知机中,我们知道一个超平面将特征空间分成两个部分,实例在不同的子空间中则被分为相对应的类。但是存在一个问题在于,我们不知道一个新输入的实例,它属于一个类的概率是多少。例如新输入的实例非常接近超平面,它被分为A类的概念为51%,分为B类的概念为49%。在感知机中将它分为了A类,但是为给出概念。
为了得到这一概率,我们引入了Sigmoid函数:
s
i
g
m
o
i
d
(
x
)
=
1
1
+
e
−
x
sigmoid(x)=\frac{1}{1+e^{-x}}
sigmoid(x)=1+e−x1
Sigmoid函数能够将线性回归产生的值
(
−
∞
,
+
∞
)
(-\infty,+\infty)
(−∞,+∞)转换到
(
0
,
1
)
(0,1)
(0,1)区间内,而概率的取值也在
(
0
,
1
)
(0,1)
(0,1)内,这样,就可以显示一个实例被分为一个类的概率是多少了。
二项逻辑斯谛回归
公式推导
逻辑斯蒂函数的一般形式,其分布具有以下分布函数和密度函数:
F
(
x
)
=
P
(
X
≤
x
)
=
1
1
+
e
−
(
x
−
u
)
/
r
F(x)=P(X\leq x) =\frac{1}{1+e^{-(x-u)/r}}
F(x)=P(X≤x)=1+e−(x−u)/r1
f ( x ) = F ′ ( x ) = 1 ( 1 + e − ( x − u ) / r ) 2 × e − ( x − u ) / r × 1 r = e − ( x − u ) / r r ( 1 + e − ( x − u ) / r ) 2 f(x)=F'(x)=\frac{1}{(1+e^{-(x-u)/r})^2}\times e^{-(x-u)/r}\times\frac{1}{r}=\frac{e^{-(x-u)/r}}{r(1+e^{-(x-u)/r})^2} f(x)=F′(x)=(1+e−(x−u)/r)21×e−(x−u)/r×r1=r(1+e−(x−u)/r)2e−(x−u)/r
式中,μ为位置参数, γ > 0 \gamma>0 γ>0为形状参数。
分布函数以(μ,1/2)为中心对称,满足:
F
(
−
x
−
μ
)
−
1
2
=
−
F
(
x
+
μ
)
+
1
2
F(-x-\mu)-\frac{1}{2}=-F(x+\mu)+\frac{1}{2}
F(−x−μ)−21=−F(x+μ)+21
形状参数
γ
\gamma
γ的值越小,分布函数曲线在中心附近增长得越快。
现在,我们让
μ
\mu
μ取0,
γ
\gamma
γ取1,即得到我们在逻辑斯蒂回归中使用的函数:
F
(
x
)
=
1
1
+
e
(
−
x
)
=
1
+
e
x
e
x
F(x)=\frac{1}{1+e^{(-x)}}=\frac{1+e^x}{e^x}
F(x)=1+e(−x)1=ex1+ex
采用上式,我们将线性回归产生的值代入到sigmoid函数之中,可得:
P
(
Y
=
1
∣
x
)
=
1
1
+
e
x
p
(
−
(
w
∗
x
)
+
b
)
P(Y=1|x)=\frac{1}{1+exp(-(w*x)+b)}
P(Y=1∣x)=1+exp(−(w∗x)+b)1
P ( Y = 0 ∣ x ) = 1 − P ( Y = 1 ∣ x ) = 1 − 1 1 + e x p ( − ( w ∗ x ) + b ) = e x p ( − ( w ∗ x ) + b ) 1 + e x p ( − ( w ∗ x ) + b ) P(Y=0|x)=1-P(Y=1|x)=1-\frac{1}{1+exp(-(w*x)+b)}=\frac{exp(-(w*x)+b)}{1+exp(-(w*x)+b)} P(Y=0∣x)=1−P(Y=1∣x)=1−1+exp(−(w∗x)+b)1=1+exp(−(w∗x)+b)exp(−(w∗x)+b)
二项逻辑斯谛回归模型是一种分类模型,由条件概率分布 P ( Y ∣ X ) P(Y|X) P(Y∣X)表示。这里,随机变量x取值为实数,随机变量 γ \gamma γ取值为1或0。
这样,我们就将范围为实数的线性回归产生的值转变为逻辑斯蒂回归中仅在 ( 0 , 1 ) (0,1) (0,1)范围之内。
逻辑斯谛回归仅对二分类的问题有效,我们可以比较 P ( Y = 1 ∣ x ) P(Y=1|x) P(Y=1∣x)和 P ( Y = 0 ∣ x ) P(Y=0|x) P(Y=0∣x)两个条件概率值的大小,将实例x分到概率较大的那一类,同时也能得知分成两种类别的可能性是多少。
逻辑斯蒂回归与几率
一个事件的几率是指该事件发生的概率与该事件不发生的概率的比值。如果事件发生的概率是p,那么该事件的几率是
p
1
−
p
\frac{p}{1-p}
1−pp ,该事件的对数几率或
l
o
g
i
t
logit
logit函数是:
l
o
g
i
t
(
p
)
=
l
o
g
p
1
−
p
logit(p)=log\frac{p}{1-p}
logit(p)=log1−pp
我们将逻辑斯蒂回归的
P
P
P代入,可得:
l
o
g
P
(
Y
=
1
∣
x
)
1
−
P
(
Y
=
1
∣
x
)
=
l
o
g
e
x
p
(
−
(
w
∗
x
)
+
b
)
1
+
e
x
p
(
−
(
w
∗
x
)
+
b
)
1
−
e
x
p
(
−
(
w
∗
x
)
+
b
)
1
+
e
x
p
(
−
(
w
∗
x
)
+
b
)
=
w
∗
x
+
b
log\frac{P(Y=1|x)}{1-P(Y=1|x)}=log\frac{\frac{exp(-(w*x)+b)}{1+exp(-(w*x)+b)}}{1-\frac{exp(-(w*x)+b)}{1+exp(-(w*x)+b)}}=w*x+b
log1−P(Y=1∣x)P(Y=1∣x)=log1−1+exp(−(w∗x)+b)exp(−(w∗x)+b)1+exp(−(w∗x)+b)exp(−(w∗x)+b)=w∗x+b
通过上式我们知道,通过几率的概念对线性函数进行转换,可以得到逻辑斯蒂回归公式。
一个直观的理解是,对于上式,分子是 y = 1 y=1 y=1的概率,而分母是 y ≠ 1 y≠1 y=1的概率,显然 w x + b wx+b wx+b越大, y = 1 y=1 y=1的概率越大,也就是实例点 x x x在 y = 1 y=1 y=1的一侧距离分离超平面越远,则 y = 1 y=1 y=1的概率越大。
模型参数估计
设:
P
(
Y
=
1
∣
x
)
=
π
(
x
)
,
P
(
Y
=
0
∣
x
)
=
1
−
π
(
x
)
P(Y=1|x)=\pi(x),P(Y=0|x)=1-\pi(x)
P(Y=1∣x)=π(x),P(Y=0∣x)=1−π(x)
似然函数为:
∏
i
=
1
N
[
π
(
x
i
)
]
y
i
[
1
−
π
(
x
i
)
]
1
−
y
i
\prod_{i=1}^N[\pi(x_i)]^{y_i}[1-\pi{(x_i)}]^{1-y_i}
i=1∏N[π(xi)]yi[1−π(xi)]1−yi
为了计算方便,我们对似然函数取对数,得到对数似然函数:
L
(
w
)
=
∑
i
=
1
N
[
y
i
l
o
g
π
(
x
i
)
+
(
i
−
y
i
)
l
o
g
(
1
−
π
(
x
i
)
)
]
=
∑
i
=
1
N
[
y
i
l
o
g
π
(
x
i
)
1
−
π
(
x
i
)
+
l
o
g
(
1
−
π
(
x
i
)
)
]
=
∑
i
=
1
N
[
y
i
(
w
∗
x
i
)
−
l
o
g
(
1
+
e
x
p
(
w
∗
x
i
)
)
]
\begin{aligned} L(w) & =\sum_{i=1}^N[y_ilog\pi (x_i)+(i-y_i)log(1-\pi(x_i))] \\ & =\sum_{i=1}^N[y_ilog\frac{\pi(x_i)}{1-\pi(x_i)}+log(1-\pi(x_i))] \\ & =\sum_{i=1}^N[y_i(w*x_i)-log(1+exp(w*x_i))] \end{aligned}
L(w)=i=1∑N[yilogπ(xi)+(i−yi)log(1−π(xi))]=i=1∑N[yilog1−π(xi)π(xi)+log(1−π(xi))]=i=1∑N[yi(w∗xi)−log(1+exp(w∗xi))]
损失函数
交叉熵的定义:
H
(
x
,
y
)
=
−
∑
i
=
1
n
x
i
l
n
y
i
4
H(x,y)=-\sum_{i=1}^nx_ilny_{i^4}
H(x,y)=−i=1∑nxilnyi4
其中x、y都是表示概率分布。假设x是正确的概率分布,而y是我们预测出来的概率分布,这个公式算出来的结果,表示y与正确答案x之间的错误程度(即:y错得有多离谱),结果值越小,表示y越准确,与x越接近。
因此逻辑斯蒂的损失函数可以看做是交叉熵损失函数。
那么为什么线性回归我们使用均方误差作为损失函数,而逻辑斯蒂回归选取交叉熵作为损失函数呢?这是因为均方误差损失函数会存在如下缺陷:
- 如果将均方误差损失函数其用于逻辑回归,则为参数的非凸函数。只有当函数为凸函数时,梯度下降才收敛到全局最小值。
- 对于均方误差损失函数来说,在更新w时,w的梯度跟激活函数的梯度成正比,激活函数梯度越大,w调整就越快,训练收敛就越快,但是Simoid函数在值非常高时候,梯度是很小的,比较平缓。而交叉熵损失函数在更新w时,w的梯度跟激活函数的梯度没有关系了,不存在梯度消失的现象。
多项逻辑斯蒂回归
多项逻辑斯蒂回归:
P
(
Y
=
k
∣
x
)
=
e
x
p
(
w
k
∗
x
)
1
+
∑
k
=
1
K
−
1
e
x
p
(
w
k
∗
x
)
,
k
=
1
,
2
,
.
.
.
,
K
−
1
P(Y=k|x)=\frac{exp(w_k*x)}{1+\sum_{k=1}^{K-1}exp(w_k*x)},k=1,2,...,K-1
P(Y=k∣x)=1+∑k=1K−1exp(wk∗x)exp(wk∗x),k=1,2,...,K−1
P ( Y = k ∣ x ) = 1 1 + ∑ k = 1 K − 1 e x p ( w k ∗ x ) P(Y=k|x)=\frac{1}{1+\sum_{k=1}^{K-1}exp(w_k*x)} P(Y=k∣x)=1+∑k=1K−1exp(wk∗x)1
参数求解的方法一样可以采用极大似然估计法。
代码
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn import datasets
# 加载鸢尾花数据集
iris = datasets.load_iris()
df = pd.DataFrame(iris.data)
df['label'] = iris.target
df.columns = ['sepal_length', 'sepal_width', 'petal_length', 'petal_width', 'label']
df.head()
# 加载鸢尾花数据集
x = iris.data
y = iris.target
# 加载鸢是三分类,将它转为二分类
y = np.array([1 if i == 1 else 0 for i in y])
# 划分训练集和测试集
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=100)
x_train,y_train
class Model:
def __init__(self):
self.w = 0 #系数
self.b = 0 #截距
self.trainSet = 0 #训练集特征
self.label = 0 #训练集标签
self.learning_rate = None #学习率
self.n_iters = None #实际迭代次数
self.accurancy = None #准确率
self.tol = 1.0e-4 #停止迭代的容忍度
self.llList = [] #记录似然值的列表
def train(self, x, y, method, n_iters=1000, learning_rate=0.01):
self.trainSet = x
self.label = y
# 梯度下降
if method.lower() == "gradient":
self.__train_gradient(n_iters, learning_rate)
elif method.lower() == "newton":
# 拟牛顿法
self.__train_newton(n_iters)
#求p(y=1|x)以及似然值LL
def PVandLLV(self, X, Y, W):
wx = np.dot(X, W.T)
p_value = np.exp(wx) / (1 + np.exp(wx))
LLarray = -1.*np.multiply(Y, wx) + np.log(1 + np.exp(wx))
return p_value, LLarray.sum()
def __calGradient(self, X, Y, Ypre):
"""calculate Gradient Matrix"""
Gw = -1.*np.multiply((Y - Ypre), X).sum(axis=0)
return Gw
def __calHessian(self, X, Ypre):
"""calculate Hessian Matrix"""
Hw = np.dot(np.dot(X.T, np.dot(np.diag(Ypre.reshape(-1)), np.diag(1-Ypre.reshape(-1)))), X)
#为了更直观的理解,展示下拆解开求解的方法
#Hw = np.zeros((X.shape[1], X.shape[1]))
#for i in range(n_samples):
# xxt = np.dot(X[i,:].reshape(-1,1),X[i,:].reshape(1,-1))
# Hw += xxt*Ypre[i]*(1-Ypre[i])
return Hw
#训练,梯度下降法
def __train_gradient(self, n_iters, learning_rate):
n_samples, n_features = self.trainSet.shape
X = self.trainSet
y = self.label
#合并w和b,在X尾部添加一列全是1的特征
X2 = np.hstack((X, np.ones((n_samples, 1))))
#将y转置变为(n_samples,1)的矩阵
Y = np.expand_dims(y, axis=1)
#初始化特征系数W
W = np.zeros((1, n_features+1))
#初始化误差,更新前后的误差之差,训练次数
Ypreprob, LL0 = self.PVandLLV(X2, Y, W)
self.llList.append(LL0)
deltaLL = np.inf
n = 0
while (n<n_iters) and (LL0>self.tol) and (abs(deltaLL)>self.tol):
#计算梯度,更新W
gra = self.__calGradient(X2, Y, Ypreprob)
W = W - learning_rate*gra/n_samples
#计算更新后的误差,并留下来
Ypreprob, LL1 = self.PVandLLV(X2, Y, W)
deltaLL = LL0 - LL1
LL0 = LL1
self.llList.append(LL0)
n += 1
self.n_iters = n
self.w = W.flatten()[:-1]
self.b = W.flatten()[-1]
Ypre = np.argmax(np.column_stack((1-Ypreprob,Ypreprob)), axis=1)
self.accurancy = sum(Ypre==y)/n_samples
print("第{}次停止迭代,似然值为{},准确率为{}".format(self.n_iters, self.llList[-1], self.accurancy))
print("w:{};\nb:{}".format(self.w, self.b))
return
#训练,牛顿法
def __train_newton(self, n_iters):
n_samples, n_features = self.trainSet.shape
X = self.trainSet
y = self.label
#合并w和b,在X尾部添加一列全是1的特征
X2 = np.hstack((X, np.ones((n_samples, 1))))
#将y转置变为(n_samples,1)的矩阵
Y = np.expand_dims(y, axis=1)
#初始化特征系数W
W = np.zeros((1, n_features+1))
#初始化误差,更新前后的误差之差,训练次数
Ypreprob, LL0 = self.PVandLLV(X2, Y, W)
self.llList.append(LL0)
deltaLL = np.inf
n = 0
while (n<n_iters) and (LL0>self.tol) and (abs(deltaLL)>self.tol):
Gw = self.__calGradient(X2, Y, Ypreprob)
Hw = self.__calHessian(X2, Ypreprob)
W = W - np.dot(Gw, np.linalg.pinv(Hw))
#计算更新后的误差,并留下来
Ypreprob, LL1 = self.PVandLLV(X2, Y, W)
deltaLL = LL0 - LL1
LL0 = LL1
self.llList.append(LL0)
n += 1
self.n_iters = n
self.w = W.flatten()[:-1]
self.b = W.flatten()[-1]
Ypre = np.argmax(np.column_stack((1-Ypreprob,Ypreprob)), axis=1)
self.accurancy = sum(Ypre==y)/n_samples
print("第{}次停止迭代,似然值为{},准确率为{}".format(self.n_iters, self.llList[-1], self.accurancy))
print("w:{};\nb:{}".format(self.w, self.b))
return
#自编的梯度下降法进行拟合
logit_gd = Model()
logit_gd.train(x, y, method="gradient", n_iters=5000, learning_rate=0.05)
plt.plot(range(logit_gd.n_iters+1), logit_gd.llList)
plt.show()