目录
这一系列是学习公众号“机器学习实验室”的笔记,跟着大佬的脚步一个个实现,逻辑回归和线性回归的实现很像,主要是公式推导那有些费时间。
数学推导
notations
X
=
{
x
i
⃗
,
i
=
1
,
2
,
.
.
.
,
m
}
X
⊆
R
m
×
d
X=\{\vec{x_i},i=1,2,...,m\} \quad X\subseteq \mathbb{R}^{m\times d}
X={xi,i=1,2,...,m}X⊆Rm×d
x
i
⃗
=
(
x
i
1
,
x
i
2
,
.
.
.
,
x
i
d
)
x
i
⊆
R
1
×
d
\quad \vec{x_i}=(x_{i1},x_{i2},...,x_{id}) \quad x_i\subseteq \mathbb{R}^{1\times d}
xi=(xi1,xi2,...,xid)xi⊆R1×d
Y
=
{
y
i
,
i
=
1
,
2
,
.
.
.
,
m
}
Y
⊆
R
m
×
1
Y=\{y_i,i=1,2,...,m\} \quad Y\subseteq \mathbb{R}^{m\times 1}
Y={yi,i=1,2,...,m}Y⊆Rm×1
y
i
⃗
⊆
R
\quad \vec{y_i}\subseteq \mathbb{R}
yi⊆R
objective function
最大化似然函数:
( W ⋆ , b ⋆ ) = a r g m a x ( w , b ) ∑ i = 1 m log p ( y i ∣ x i ; w , b ) (W^{\star},b^{\star})=\mathop{argmax}\limits_{(w,b)}\sum_{i=1}^{m}\log{p(y_i|x_i;w,b)} (W⋆,b⋆)=(w,b)argmax∑i=1mlogp(yi∣xi;w,b)
思想推导
广义线性模型:
y = g − 1 ( w T x + b ) y=g^{-1}(w^Tx+b) y=g−1(wTx+b),这个函数是可微的,有 f = g − 1 f=g^{-1} f=g−1。
我们取sigmoid函数
y
=
1
1
+
e
−
x
y=\frac{1}{1+e^{-x}}
y=1+e−x1
注意:对数线性回归不等于对数几率回归。
对数线性回归 | 对数几率回归 |
---|---|
log-linear regression | logistic regression |
l n y = w T x + b ln{y}=w^Tx+b lny=wTx+b | l n y 1 − y = w T x + b ln{\frac{y}{1-y}=w^Tx+b} ln1−yy=wTx+b |
y = e w T x + b y=e^{w^Tx+b} y=ewTx+b | y = 1 1 + e − ( w T x + b ) y=\frac{1}{1+e^{-(w^{T}x+b)}} y=1+e−(wTx+b)1 |
推导还是比较简单的:
y
=
1
1
+
e
−
(
w
T
x
+
b
)
y=\frac{1}{1+e^{-(w^{T}x+b)}}
y=1+e−(wTx+b)1
1 + e − ( w T x + b ) = 1 y 1+e^{-(w^Tx+b)}=\frac{1}{y} 1+e−(wTx+b)=y1
e − ( w T x + b ) = 1 − y y e^{-(w^Tx+b)}=\frac{1-y}{y} e−(wTx+b)=y1−y
两边去对数:
w T x + b = l n y 1 − y w^Tx+b=ln{\frac{y}{1-y}} wTx+b=ln1−yy
记 y y y 为类后验概率 p ( y = 1 ∣ x ) p(y=1|x) p(y=1∣x), 1 − y 1-y 1−y 则为 p ( y = 0 ∣ x ) p(y=0|x) p(y=0∣x)。
-
对一个样本来说, p ( y i ∣ x i ) = y i ∗ p ( y = 1 ∣ x i ) + ( 1 − y i ) ∗ p ( y = 0 ∣ x i ) p(y_i|x_i)=y_i*p(y=1|x_i)+(1-y_i)*p(y=0|x_i) p(yi∣xi)=yi∗p(y=1∣xi)+(1−yi)∗p(y=0∣xi),标签为1时,计算 p ( y = 1 ∣ x ) p(y=1|x) p(y=1∣x),标签为0时,计算 p ( y = 0 ∣ x ) p(y=0|x) p(y=0∣x)。
-
对于整个样本集来说,想要求使得样本联合概率最大的参数。 ( W ⋆ , b ⋆ ) = a r g m a x ( w , b ) L ( w , b ) (W^{\star},b^{\star})=\mathop{argmax}\limits_{(w,b)}L(w,b) (W⋆,b⋆)=(w,b)argmaxL(w,b)
L ( w , b ) = ∏ i = 1 m p w , b ( y i ∣ x i ) L(w,b)=\prod_{i=1}^{m}p_{w,b}(y_i|x_i) L(w,b)=∏i=1mpw,b(yi∣xi)
事情就变得很奇妙了,取对数后有
l n L ( w , b ) = ∑ i = 1 m ln p w , b ( y i ∣ x i ) ln{L(w,b)}=\sum_{i=1}^{m}\ln{p_{w,b}(y_i|x_i)} lnL(w,b)=∑i=1mlnpw,b(yi∣xi)
看一下log函数,发现是单调递增的,也就是说最大化联合概率就转换成最大化每个样本概率的对数和。
这就是最大似然概率的公式,最大似然法的优化思想。
3. 把1带入2的
L
(
w
,
b
)
L(w,b)
L(w,b) 有
∑
i
=
1
m
(
y
i
(
w
T
x
+
b
)
−
l
n
(
1
+
e
w
T
x
+
b
)
)
\sum_{i=1}^{m}{(y_i(w^Tx+b)-ln(1+e^{w^Tx+b}))}
∑i=1m(yi(wTx+b)−ln(1+ewTx+b)),最大化这个式子就等于最小化
∑ i = 1 m ( − y i ( w T x + b ) l n ( 1 + e w T x + b ) ) \sum_{i=1}^{m}{(-y_i(w^Tx+b)ln(1+e^{w^Tx+b}))} ∑i=1m(−yi(wTx+b)ln(1+ewTx+b))。
partial derivative
∂ E ( w , b ) ∂ w = ∑ i = 1 m − y i x i + x i e w T x + b 1 + e w T x + b \frac{\partial{E(w,b)}}{\partial{w}}=\sum_{i=1}^{m}{-y_ix_i+x_i\frac{e^{w^Tx+b}}{1+e^{w^Tx+b}}} ∂w∂E(w,b)=∑i=1m−yixi+xi1+ewTx+bewTx+b
∂ E ( w , b ) ∂ b = ∑ i = 1 m − y i + e w T x + b 1 + e w T x + b \frac{\partial{E(w,b)}}{\partial{b}}=\sum_{i=1}^{m}-y_i+\frac{e^{w^Tx+b}}{1+e^{w^Tx+b}} ∂b∂E(w,b)=∑i=1m−yi+1+ewTx+bewTx+b
matrix version
我们利用Numpy实现,数据存放为二维数组,在代码实现的时候需要看一下矩阵相乘合不合法。
手动实现逻辑回归
首先分析有几个部分:LogisticRegression作为一个类,应该满足喂入训练数据x和y,然后拟合好一个线性函数F,参数为W和b;此外根据这个线性函数通过sigmoid可以得到测试数据属于类别的概率。
由此我们可以很快想到造一个训练器train、一个预测器predict和一个评价器evaluation。
因为思路和线性回归是一样的,仅仅改变了一下predict和evaluation函数,因此这里重点放这两个。
训练器 train
思路
- 训练器的input是x_train和y_train,output是weight和bias。
- 优化思想是最大似然法,出发点是最大化似然函数,这样loss就可以确定为
∑
i
=
1
m
log
p
(
y
i
∣
x
i
;
w
,
b
)
\sum_{i=1}^{m}\log{p(y_i|x_i;w,b)}
∑i=1mlogp(yi∣xi;w,b),推导见上。
其中 p ( y i ∣ x i ; w , b ) = y i p 1 ( x i ^ ; w , b ) + ( 1 − y i ) p 0 ( x i ^ ; w , b ) p(y_i|x_i;w,b)=y_ip_1(\hat{x_i};w,b)+(1-y_i)p_0(\hat{x_i};w,b) p(yi∣xi;w,b)=yip1(xi^;w,b)+(1−yi)p0(xi^;w,b) p 0 ( x i ^ ; w , b ) = 1 − p 1 ( x i ^ ; w , b ) p_0(\hat{x_i};w,b)=1-p_1(\hat{x_i};w,b) p0(xi^;w,b)=1−p1(xi^;w,b) - 考虑到特征维数较高,就不用正规方程解法,直接用数值优化方法——梯度下降,对loss函数求导见上,这和线性回归是一样的。
因此我们可以写出伪代码:
initialize parameters w, b
for epoch in epochs:
caculate loss, dw, db
update w with dw
update b with db
return w, b
第一步是初始化参数,权重与特征维数一致,偏差是个常数。
第二步是求loss和偏导,我们可以通过一个函数完成。这个函数要求输入x_train,y_train和参数,返回loss,dw,db。
第三步就是利用返回的偏导更新参数,从稍后的图可以看到loss会震荡,然后平稳。
代码实现
重点讲loss这块,其他都和线性回归差不多,这里好像写错了一点,我之后再修正。
# 求loss,dw,db
def loss(self, X, y, W, b):
num_train = X.shape[0]
num_feature = X.shape[1]
beta_x = np.dot(X, W) + b
y_hat = self.sigmoid(beta_x)
p_y = y * y_hat + (1 - y) * (1 - y_hat)
# loss = np.sum(( y * np.log(beta_x) - np.log(1 + beta_x)))
loss = np.sum(- y + self.sigmoid(beta_x))
dW = np.dot(X.T, - y + self.sigmoid(beta_x)) / num_feature
db = np.sum(- y + self.sigmoid(beta_x)) / num_train
# print(dW)
return y_hat, loss, dW, db
预测器 predict
思路
input是存放的使loss最小的parameters和x_test,output是计算得到的y_predict。
概率大于等于0.5判为1,否则0,直接上代码。
代码
def predict(self, X, params):
W = params['W']
b = params['b']
y_predict = self.sigmoid(np.dot(X, W) + b)
y_predict[y_predict >= 0.5] = 1
y_predict[y_predict < 0.5] = 0
评价器 evaluation
思路
input是预测得到的值y_predict和原值y_test,output是我们定义的准确率,精确率,召回率。
代码
def evalution(self, y_predict, y):
y_test = y_test.ravel()
y_predict = y_predict.ravel()
num_test = y_test.shape[0]
# succeed: 0 fail: 1,-1
outcome = y_predict - y_test
# succeed: 0 fail: 1
outcome[outcome == -1] = 1
Counter(outcome)
Counter(y_predict)
num_tp_tn = sum(outcome==0)
num_fp = np.dot(y_predict.T, outcome)
num_tp = sum(y_predict==1) - num_fp
num_tn = num_tp_tn - num_tp
num_fn = sum(y_predict==0) - num_tn
# print(outcome, y_predict, y_test)
# print(num_tp_tn, num_fp, num_tp, num_tn, num_fn)
# print(list(y_test).count(0), list(y_test).count(1))
accuracy = num_tp_tn / num_test
precision = num_tp / (num_tp + num_fp)
recall = num_tp / (num_tp + num_fn)
return accuracy, precision, recall
与sklearn比较
我加入sklearn的包来比较模型的效果,在学习率为1,迭代10000次的时候在这个测试集上达到和sklearn一样的效果。
ours: 0.9122807017543859 0.9210526315789473 0.9459459459459459
sklearn: 0.9122807017543859 0.9210526315789473 0.9459459459459459
学习率
在0.5学习率、迭代10000次的时候并不好。
学习率为1时可以看到和sklearn的预测完全一致。
损失
学习过程中loss一直在震荡,这里其实我没有完全按照loss函数(因为值溢出了),而是定义了一个Loss。
降维可视化
这是比较尴尬的,因为维数比较高,可视化的话最好是降维处理,比如主成分分析,或者用可视化包,但是我太懒了,选了第二维和第三维特征,看不出来,但我想表达如果学得好明显就能把数据在二维平面分开。
完整代码
https://github.com/oneoyz/ml_algorithms/tree/master/logistic_regression