JSong @2016.06.13
本系列文章不适合入门，是作者综合各方资源和个人理解而得. 另外最好有数学基础, 因为数学人一言不合就会上公式.
简单模型的魅力在于它能从各个角度去欣赏. 逻辑回归是最简单的二分类模型之一，实际应用中二分类最常见，如判定是否是垃圾邮件，是否是人脸，是否值得借贷等, 而概率模型对于这类问题有得天独厚的优势. 本文将从各个角度来理解逻辑回归，并指出它是一个概率模型、对数线性模型、交叉熵模型.
我们先从线性回归谈起。 考察 m 个变量和 y 之间的线性关系：
\[y\sim x_1+x_2+\cdots +x_m\]
根据要求我们需要找到 m+1 个回归系数 θ_i，使得
\[\min ||y-\hat{y}||_{F},\quad \hat{y}=h_{\theta}(x)=\sum_{i=0}^m \theta_i x_i=\theta^T x\]
其中 x_0=1.
1. 线性代数下的回归模型
如果每个变量都有 n 个样本，即 x_i ∈ R^n ，则上述问题等价于求解线性方程组：
\[X\theta =y, \quad X=(x_0,x_1,\ldots,x_m)\]
一般来讲，上述方程有唯一解(最小二乘法)：
\[\theta=X^{+}y=(X^TX)^{-1}X^Ty\]
讲到这，车就可以开始上路了。
2. 二分类问题简介
分类模型比较常见，例如是否是垃圾邮件等等。而最简单的分类就是二分类，且多分类一定程度也可以转化成二分类问题。 给定一组有标签的训练数据
\[(X,y)=\{(x^{i}\in \mathbb{R}^n, y^{i}\in \{0,1\}): i=1,2,\ldots,m\}\]
在下图中我们随机给了一组散点图，红色代表y=1，蓝色代表y=0.
我们的目的就是找到一种模型将这两类样本点分开(当然实际应用中可能没有这么好的情况。 这里只是给一个Demo)。 最直接的想法就是找一个超平面
\[\theta^T x=0\]
使得红色点和蓝色点分别分布于超平面的两端. 此时对于任意一个样本点 x^{(i)}, 不失一般性
\[\theta^T x^{(i)}\]
可以代表样本点到该超平面的距离。 因为红色样本点在超平面上方，即红点对应距离为正数，同理蓝点对应的距离为负数。 为避免区分红点或者蓝点，我们可以用
\[(2y-1)\theta^T x\]
来代替， 此时只要被正确分类， 该距离则为正数，否则为负数。 于是可得相应的损失函数为
\[\max(0,(2y-1)\theta^T x)\]
综合一下可得最简单的线性分类模型为：
\[\arg_\theta \min \sum_{i=1}^{m} \max(0,(2y_i-1)\theta^T x^{(i)})\]
这个损失函数有很多不好的地方，而且也不可导。接下来我们介绍更好的逻辑回归模型，其背后有很多的解释。
2.1 逻辑回归模型
令
\[\theta^T x=\theta_0 +\sum_{i=1}^{m} \theta_i x_i\]
且我们假定
\[P(y=1|x; \theta)=h_{\theta}(x)=g(\theta^T x)=\frac{1}{1+e^{-\theta^T x}}\]
\[P(y=0|x; \theta)=1-h_{\theta}(x)\]
其中
\[g(z)=\frac{1}{1+e^{-z}}\]
是 logistic 分布函数，它的图像和密度函数
\[e^{-z}/(1+e^{-z})^2\]
图像见下图. 原则上任何连续的分布函数都可以，其把实数映射到0到1之间. 至于为啥选择logistic函数, 我们之后再讨论.
这里我们采用极大似然估计. 由于y取值的特殊性, 上面的模型假设等价于
\[P(y|x;\theta)=(h_{\theta}(x))^y \,(1-h_{\theta}(x))^{1-y}\]
这样我们便有似然函数：
\[l(\theta)=\log L(\theta)=\log \prod_{i=1}^{m} p(y^{(i)}|x^{(i)};\theta)\]
\[=\sum_{i=1}^m [y^{(i)} \log h_{\theta}(x^{(i)})+ (1-y^{(i)}) \log (1-h_{\theta}(x^{(i)})]\]
于是得到逻辑回归模型为
\[\arg_{\theta} \max l(\theta)\]
2.2 几何意义（超平面分类）
重新考虑模型假设, 我们可以得到
\[\theta^T x=\log \frac{p}{1-p}, \quad p=P(y=1|x)\]
同样，逻辑回归模型也是用超平面来分类的. 且相应的极大似然估计可以写成
\[\min -l(\theta)= \sum_{i=1}^{m} \log[(2y-1)(y-h_{\theta}(x))]=\sum_{i=1}^{m} \log (1+e^{(2y-1)\theta^T x})\]
也即相应的损失函数为
\[\log(1+\exp((2y-1)\theta^T x))\]
图像见下方. 由统计假设可知，给定一个样本点，如果它被正确分类，则
\[(2y-1)(y-h_{\theta}(x))\]
应该越小越好, 如y=1时,h_θ(x)越接近1越好. 相应的 (2y−1)θ^Tx 也是越小越好. 又根据第一节的讨论我们有, 当某样本点被正确分类时
\[(2y-1)\theta^T x= \mbox{distance of sample point and hyperplane}\,\, \theta^T x=0 \]
此处
\[\log(1+\exp(x)) \quad and \quad \max(0,x)\]
的作用类似，但连续性更好. 另外有些文章会写逻辑回归的损失函数为
\[\log(1+\exp(-y\cdot \theta^T x))\]
这是因为该损失函数中的 y 取值范围为-1和1.
2.3 线性代数意义
在上一节，我们给出了表达式
\[\theta^T x=\log \frac{p}{1-p}, \quad p=P(y=1|x)\]
可以看出它相当于把 y 映射到了[0,1]区间。而且 p/(1-p) 是事件发生与事件不发生的概率之比，称为事件的发生比 (the odds of experiencing an event), 简称为 odds 。
2.4 信息论解释(交叉熵模型)
令随机变量
\[p\in \{y,1-y\}, \quad q \in\{\hat{y},1-\hat{y}\}\]
其中 \hat{y} 为模型拟合得到的y, 则它们之间的交叉熵 (cross entropy) 为
\[H(p,q)= -\sum_i p_i \log q_i =-[y\log \hat{y}+(1-y)\log(1-\hat{y})]\]
我们知道熵常用于度量一个随机变量所包含的信息量， 而交叉熵可以用来度量两个随机变量之间的相似性, 从上面式子可以看出逻辑回归的极大似然等价于最小化交叉熵.
事实上真正等价的是交叉熵与极大似然估计，有兴趣的同学可以自己证明。
2.5 神经网络模型
这个就不多说了，逻辑回归是一个最简单的神经网络模型
2.6 梯度下降法参数求解
按照极大似然的来求导：
\[\bigtriangledown_{\theta} l(\theta)=\sum_{i=1}^{m} (y^{(i)}-h_{\theta}(x^{(i)}))x^{(i)}= X^T[Y-h_{\theta}(X)]\]
其中 x^{(i)} 在 X 中作为行向量存储. 于是可得参数 θ 的迭代式：
\[ \theta:=\theta+\alpha X^T[Y-h_{\theta}(X)] \]
其中 α 为步长, 另外因为是最大化 l(θ) ，所以是沿着梯度的正方向寻找。
1.7 python代码
import numpy as np
import matplotlib.pyplot as plt
import time
%pylab inline# calculate the sigmoid function
def sigmoid(inX):return 1.0 / (1 + np.exp(-inX))def trainLogRegres(train_x, train_y, opts):# train a logistic regression model using some optional optimize algorithm# train_x is a mat datatype, each row stands for one sample# train_y is mat datatype too, each row is the corresponding label# opts is optimize option include step and maximum number of iterations# calculate training timestartTime = time.time()train_x=np.asmatrix(train_x)train_y=np.asmatrix(train_y)numSamples, numFeatures =train_x.shapealpha = opts['alpha']; maxIter = opts['maxIter']weights = np.ones((numFeatures, 1))for k in range(maxIter):err = train_y - sigmoid(train_x * weights)weights = weights + alpha * train_x.T * errprint 'Congratulations, training complete! Took %fs!' % (time.time() - startTime)return weights# show your trained logistic regression model only available with 2-D data
def showLogRegres(weights, train_x, train_y):# notice: train_x and train_y is mat datatypenumSamples, numFeatures = train_x.shapeif numFeatures != 3:print "Sorry! I can not draw because the dimension of your data is not 2!"return 1# draw all samples idx1=np.asarray(train_y)==1idx2=np.asarray(train_y)==0plt.plot(x[:,1][idx1].T,x[:,2][idx1].T,'or')plt.plot(x[:,1][idx2].T,x[:,2][idx2].T,'ob')# draw the classify linemin_x = min(train_x[:, 1])[0, 0]max_x = max(train_x[:, 1])[0, 0]weights = weights.getA() # convert mat to arrayy_min_x = float(-weights[0] - weights[1] * min_x) / weights[2]y_max_x = float(-weights[0] - weights[1] * max_x) / weights[2]plt.plot([min_x, max_x], [y_min_x, y_max_x], '-g')plt.xlabel('X1'); plt.ylabel('X2')# 测试样例
len_samples=500
x=np.random.randn(len_samples,2)
idx1=x[:,0]+x[:,1]>0.5
y=np.zeros((len_samples,1))
y[idx1,0]=1
opt={'alpha':0.1,'maxIter':1000}
x=np.hstack((np.ones((len_samples,1)),x))
x=np.asmatrix(x)
y=np.asmatrix(y)
w=trainLogRegres(x,y,opt)
print w
showLogRegres(w,x,y)输出如下：
Populating the interactive namespace from numpy and matplotlib
Congratulations, training complete! Took 0.047000s!
[[-17.37216188][ 34.35470115][ 34.98381818]]
这只是最简单的样例，在实际应用中我们还可以添加惩罚项
\[ \min_{w,c} \|w\|_2+ C\sum_{i=1}^{m} \log(1+\exp(-y^{(i)}(w\cdot x^{(i)}+c)))\]
或者 L_1 范数
\[ \min_{w,c} \|w\|_1+ C\sum_{i=1}^{m} \log(1+\exp(-y^{(i)}(w\cdot x^{(i)}+c)))\]
在样本很稀疏的时候，惩罚项很有用. 另外当样本量很大，我们又极其要求速度的时候，梯度下降法也要改进，换成SGD（stochastic gradient descent ）、拟牛顿法、AGD等等
专注写文二十年，原文首发于公众号，欢迎关注
