Logistic Regression 原理及推导 python实现

一、问题引入

首先，Logistic回归是一种广义的线性回归模型，主要用于解决二分类问题。比如，现在我们有N个样本点，每个样本点有两维特征x1和x2，在直角坐标系中画出这N个样本的散点图如下图所示，



蓝色和红色分别代表两类样本。现在我们的目标是，根据这N个样本所表现出的特征以及他们各自对应的标签，拟合出一条直线对两类样本进行分类，直线的上侧属于第一类，直线的下侧属于第二类。那么我们如何寻找这条直线呢？我们知道，二维空间的一条直线可以用简单的公式表示

y=b+θ1x1+θ2x2=θTx+b

参数θT和b的选择决定了直线的位置，如果我们选取了一组参数θ和b导致直线的位置是这样的



那肯定不合理，因为两类样本完全没有被分开，而如果我们得到了这样一条直线



两类样本虽然存在一些小错误，但是基本上被分开了，所以决定这条直线位置的参数和θ 和b就是我们想要的答案。由此，我们可以看到，Logistic Regression问题最终变成了求解参数θ的问题。

二、原理

总的来说，Logistic回归的过程就是根据样本求解分类直线（或者超平面）的参数的过程，求解的方法是极大似然估计法，但是因为似然方程最后求导无法得到解析解，所以用了梯度下降法去逐步优化得到极值。

为什么要用极大似然估计法来求解参数呢？

首先，假设样本有k维特征，x={x1,x2,...,xk}, y∈{0,1}，用极大似然估计就是在求解怎样的θ 和b可以使得出现了某个x时，出现对应的y的概率最大。

然后，假设这个概率服从的是Sigmoid分布（如下图所示）。样本的每一维特征的取值在经过参数θ线性组合之后取值范围是实数集（-inf, inf），通过Sigmoid变换取值变成了（0，1）。

Sigmoid函数：

g(z)=11+e−z



三、极大似然估计求解迭代方程

首先，把问题变成一个概率问题：

在某个x和θ的取值下，y=1的概率为hθ(x)，P(y=1|x;θ)=hθ(x)

在某个x和θ的取值下，y=0的概率为1−hθ(x)，P(y=0|x;θ)=1−hθ(x)

由于y只有两种取值：0和1，因此综合两种情况，对于每一个样本点来说，P(y|x;θ)=(hθ(x))y(1−hθ(x))1−y

考虑样本集中的所有样本点，由于每个样本之间相互独立，因此它们的联合分布等于各自边际分布之积，

L(θ)=∐i=1mP(yi|xi;θ)=(hθ(xi))yi(1−hθ(xi))1−yi

这就是我们求解θ需要的似然函数，我们通过他来求解在θ为何值时，x取某个值出现某个y的概率最大。

对L(θ)取对数，因为ln(x)和x单调性相同

l(θ)=lnL(θ)=∑i=1m(yilnhθ(xi)+(1−yi)ln(1−hθ(xi)))

给出损失函数J(θ)=−1ml(θ)，对J(θ)求偏导，



理应令求偏导后的表达式等于零求极值，但是无法解析求解，因此用梯度下降法逐渐迭代，找到局部最优解。

四、梯度下降法求解局部极值

为什么梯度下降法能够做到呢？



可以看到θ的取值和J(θ)存在着一一对应的关系，让θ沿着J(θ)梯度的方向减小，可以最快速的逼近J(θ)的最小值，但其实往往找到的是极小值，局部最优。



由此可以得到θ的更新公式



五、Logistic回归算法步骤

初始化回归系数θ为1

重复下面步骤直至收敛

{

计算整个数据集的梯度

使用α x gradient更新回归系数θ

}

返回回归系数θ

六、python实现

现在，我们来解决一个实际问题：分类sklearn的make_moon数据集。该数据集每个样本有两维特征，我们选取前250个样本，绘制Logistic回归的训练过程。

#!/bin/bash/env python
#-*- coding: utf-8
__author__ = 'ZhangJin'

import numpy as np
from sklearn.datasets import make_moons

class LR:
def __init__(self):
self.dim = 2
#       self.w = np.random.random(self.dim)
self.w = np.array([1.0,1.0])
self.b = 0
self.eta = 0.2

def sigmoid(self, x):
return 1.0/(1+np.exp(-x))

def logistic_regression(self,x,y,eta):
itr = 0
self.eta = eta
row, column = np.shape(x)
xpts = np.linspace(-1.5, 2.5)
while itr <= 100:
fx = np.dot(self.w, x.T) + self.b
hx = self.sigmoid(fx)
t = (hx-y)
s = [[i[0]*i[1][0],i[0]*i[1][1]] for i in zip(t,x)]
gradient_w = np.sum(s, 0)/row * self.eta
gradient_b = np.sum(t, 0)/row * self.eta
self.w -= gradient_w
self.b -= gradient_b
ypts = (lr.w[0] * xpts + lr.b) / (-lr.w[1])
if itr%20 == 0:
plt.figure()
for i in range(250):
plt.plot(x[i, 0], x[i, 1], col[y[i]] + 'o')
plt.ylim([-1.5,1.5])
plt.plot(xpts,ypts, 'g-', lw = 2)
plt.title('eta = %s, Iteration = %s\n' % (str(eta), str(itr)))
plt.savefig('p_N%s_it%s' % (str(row), str(itr)), dpi=200, bbox_inches='tight')
itr += 1

if __name__ == '__main__':
import matplotlib.pyplot as plt
x, y = make_moons(250, noise=0.25)
col = {0:'r',1:'b'}
lr = LR()
lr.logistic_regression(x,y,eta=1.2)

得到的结果，



七、数据线性不可分的思考

这个例子中用到的make_moons数据集其实是线性不可分的。那怎么用logistic回归处理这样的线性不可分问题呢？

可以考虑把数据变换到高维（多项式变换），让数据变得可分。这里我还没有学习SVM的核技巧，后面学到了再加进我的思考。

目前每个样本只有两维x1,x2，变换后得到比如五维x1,x2,x1x2,x21,x22。升维之后用分类算法训练模型的时候参数的数量会很大，计算大量消耗资源。可以考虑用以下步骤减少样本的数量：

1. 用K-means聚若干个类簇

2. 对每个中心点，如果它周围距离小于ϵ 的点中，两类样本的数目大于某个阈值（比如2:3），就认为这个中心点处在决策边界上

3. 挑选出在决策边界上的中心点作为新的样本集，对其进行升维，分类

class kmeans:
def __init__(self, k, N, x):
self.k = k
self.N = N
# randomly choose k centroids in N points
self.centroids = x[random.sample(range(N),k),:]

def kma(self, x):
# k-means algorithm
itr = 0
while itr <= 200:
cluster_label = np.zeros([N,1])
for i in range(self.N):
temp = np.dot(np.ones((self.k,1)), np.reshape(x[i],(1,2)))
dist = np.sum((temp - self.centroids) ** 2, 1)
cluster_label[i] = np.argmin(dist)
# 更新质心
for i in range(self.k):
t = x[np.where(cluster_label==i),:][0]
self.centroids[i] = np.sum(t, 0)/np.shape(t)[0]
itr += 1

# 用K-Means找到中心点之后，需要找到边界处的中心点。
# 对于每一个中心点，如果周围M个点中两类点的比例大于3／7，那么这个中心点就是边界处的中心点
def find_border_centroids(self, epsilon = 1):
self.border_centroids = []
for i in range(self.k):
temp = np.dot(np.ones((self.N,1)), np.reshape(self.centroids[i],(1,2)))
dist = np.sum((x-temp)**2,1)
k = np.histogram(y[np.where(dist < epsilon)],bins=2)[0]
if min(k)/float(max(k)) >= 0.43 :
self.border_centroids.append(self.centroids[i])

def plot_centroids(self):
#       for i in range(N):
#           plt.plot(x[i, 0], x[i, 1], col[y[i]] + 'o')

n = (np.shape(self.border_centroids))[0]
for i in range(n):
plt.plot(self.border_centroids[i][0], self.border_centroids[i][1], 'r+')
plt.show()

此时数据集的点有1000个



1000个点用K-means找到100个聚类中心



找到决策边界的中心点



（代码可以点击此处下载）(https://github.com/zjsghww/MachineLearning)