【机器学习系列2】逻辑回归LogisticRegression——数学推导和纯Numpy实现

这一系列是学习公众号“机器学习实验室”的笔记，跟着大佬的脚步一个个实现，逻辑回归和线性回归的实现很像，主要是公式推导那有些费时间。

数学推导

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
$\overrightarrow{x_i}=(x_{i1},x_{i2},...,x_{id})x_i\subseteq {\mathbb{R}}^{1\times 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
$\overrightarrow{y_i}\subseteq \mathbb{R}$

objective function

最大化似然函数：

$(W^{\star },b^{\star })=\underset{(w,b)}{\mathrm{argmax}}{\sum }_{i=1}^m\log p(y_i|x_i;w,b)$

思想推导

广义线性模型：

$y=g^{-1}(w^{T}x+b)$，这个函数是可微的，有$f=g^{-1}$。

我们取sigmoid函数 $y=\frac{1}{1+e^{-x}}$

注意：对数线性回归不等于对数几率回归。

对数线性回归	对数几率回归
log-linear regression	logistic regression
$lny=w^{T}x+b$	$ln\frac{y}{1-y}=w^{T}x+b$
$y=e^{w^{T}x+b}$	$y=\frac{1}{1+e^{-(w^{T}x+b)}}$

推导还是比较简单的：
$y=\frac{1}{1+e^{-(w^{T}x+b)}}$

$1+e^{-(w^{T}x+b)}=\frac{1}{y}$

$e^{-(w^{T}x+b)}=\frac{1-y}{y}$

两边去对数：

$w^{T}x+b=ln\frac{y}{1-y}$

记 $y$ 为类后验概率 $p(y=1|x)$，$1-y$ 则为 $p(y=0|x)$。

1. 对一个样本来说，$p(y_i|x_i)=y_i\ast p(y=1|x_i)+(1-y_i)\ast p(y=0|x_i)$，标签为1时，计算$p(y=1|x)$，标签为0时，计算$p(y=0|x)$。

2. 对于整个样本集来说，想要求使得样本联合概率最大的参数。$(W^{\star },b^{\star })=\underset{(w,b)}{\mathrm{argmax}}L(w,b)$

   $L(w,b)={\prod }_{i=1}^m{p}_{w,b}(y_i|x_i)$

   事情就变得很奇妙了，取对数后有

   $lnL(w,b)={\sum }_{i=1}^m\ln p_{w,b}(y_i|x_i)$

   看一下log函数，发现是单调递增的，也就是说最大化联合概率就转换成最大化每个样本概率的对数和。

   这就是最大似然概率的公式，最大似然法的优化思想。

3. 把1带入2的 $L(w,b)$ 有 ${\sum }_{i=1}^m(y_i(w^{T}x+b)-ln(1+e^{w^{T}x+b}))$，最大化这个式子就等于最小化

${\sum }_{i=1}^m(-y_i(w^{T}x+b)ln(1+e^{w^{T}x+b}))$。

partial derivative

$\frac{\partial E(w,b)}{\partial w}={\sum }_{i=1}^m-y_i{x}_i+x_i\frac{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^{T}x+b}}{1+e^{w^{T}x+b}}$

matrix version

我们利用Numpy实现，数据存放为二维数组，在代码实现的时候需要看一下矩阵相乘合不合法。

手动实现逻辑回归

首先分析有几个部分：LogisticRegression作为一个类，应该满足喂入训练数据x和y，然后拟合好一个线性函数F，参数为W和b；此外根据这个线性函数通过sigmoid可以得到测试数据属于类别的概率。
由此我们可以很快想到造一个训练器train、一个预测器predict和一个评价器evaluation。
因为思路和线性回归是一样的，仅仅改变了一下predict和evaluation函数，因此这里重点放这两个。

训练器 train

思路

1. 训练器的input是x_train和y_train，output是weight和bias。
2. 优化思想是最大似然法，出发点是最大化似然函数，这样loss就可以确定为${\sum }_{i=1}^m\log p(y_i|x_i;w,b)$，推导见上。
   其中$p(y_i|x_i;w,b)=y_i{p}_1(\widehat{x_i};w,b)+(1-y_i)p_0(\widehat{x_i};w,b)$$p_0(\widehat{x_i};w,b)=1-p_1(\widehat{x_i};w,b)$
3. 考虑到特征维数较高，就不用正规方程解法，直接用数值优化方法——梯度下降，对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