逻辑回归(Logistic Regression)

逻辑回归(Logistic Regression)

1 极大似然估计(maximum likelihood estimation)

概念: 极大似然估计是一种概率论在统计学的应用，是参数评估的方法之一。

假设 已知某个样本满足满足某种概率分布，但是其中具体的参数并不清楚，参数估计通过若干次试验，观察其结果，利用结果推出参数的最大概率值。

极大似然估计就是建立在这样的思想上的：已知某个参数能使这个样本出现的概率最大(这句话很重要，使得样本同时出现的概率最大，表现为:x,y同时出现的概率最大)。如此，不会再去选择其他概率小的样本，所以把这个参数作为估计的真实值。
求最大似然函数估计值的一般步骤：

写出似然函数；
对似然函数取对数，并整理；
求导数，令导数为0，得到似然方程；
解似然方程，得到的参数即为所求；

2 Logistic算法

Logistic算法是Softmax算法的特例，按照由特殊到一般的原则，先讲逻辑斯蒂回归，再推出Softmax

2.1 算法的优缺点

优点：计算代价不高，易于理解和实现，
缺点：容易欠拟合，分类精度可能不高
使用数据类型，数值型和标称型数据

2.2 算法理解&推导

2.2.1 Sigmoid函数

Assume：
问题：
希望有这么一个f(x)，能够具有很好的逻辑判断性质，最好是能够直接表达具有特征x的样本被分到某类的概率。比如f(x)>0.5的时候能够表示x被分为正类，f(x)<0.5表示分为反类。而且我们希望f(x)总在[0, 1]之间。有这样的函数吗？
先来直观了解一下sigmoid函数：

sigmoid(x)=11+e−x $$sigmoid(x)=\frac{1}{1+e^{-x}}$$

由sigmoid的函数可得：sigmoid定义域为全体实数，值域在(0, 1)之间，并且当 x=0 $x=0$ 时 sigmoid(0)=0.5 $sigmoid(0)=0.5$ 。sigmoid看起来很像一个 阶跃函数（阶跃函数是奇异函数，t<0时，函数值为0；t=0时，函数值为1/2,；t>0时，函数值为1），但其实也不太像。

2.3 Logistic算法的推导

假设我们的样本是{ x , y}，y为1或者0，表示正类或者负类， x 是m维的样本特征向量。那么这个样本 x 属于正类，也就是y=1的“概率”可以通过下面的逻辑函数来表示：

p(y=1|x;θ)=δ(θTx)=11+e−θTx=eθTx1+eθTx $$p(y=1|x;\theta )=\delta ({\theta }^{T}x)=\frac{1}{1+e^{-{\theta }^{T}x}}=\frac{e^{{\theta }^{T}x}}{1+e^{{\theta }^{T}x}}$$

p(y=0|x;θ)=1−p(y=1|x;θ)=11+eθTx $$p(y=0|x;\theta )=1-p(y=1|x;\theta )=\frac{1}{1+e^{{\theta }^{T}x}}$$

θ 是模型参数，也就是回归系数，σ是sigmoid函数。实际上这个函数是由下面的对数几率（也就是 x 属于正类的可能性和负类的可能性的比值的对数）变换得到的：

logit(x)=ln(p(y=1|x;θp(y=0|x;θ)) $$\log it(x)=\ln (\frac{p(y=1|x;\theta }{p(y=0|x;\theta )})$$

---

举个例子(网上看到的，觉得比较好就直接拿来了)：
假设：y是关系的变量，x是影响y或者决定y的因素。例如一个女孩子喜不喜欢你，与多个自变量（因素）有关，如你人品怎样、车子是两个轮的还是四个轮的、长得胜过潘安还是形如犀利哥、有千尺豪宅还是三寸茅庐等等，我们把这些因素表示为 x1,x2,…,xm $x_1,x_2,\dots ,x_m$。
问题来了：
那这个女生怎样考量这些因素呢？
最快的方式就是把这些因素的得分都加起来，最后得到的和越大，就表示越喜欢。
但每个人心里其实都有一杆称，每个人考虑的因素不同，萝卜青菜，各有所爱嘛。例如这个女生更看中你的人品，人品的权值是0.6，不看重你有没有钱，没钱了一起努力奋斗，那么有没有钱的权值是0.001等等。我们将这些对应 x1,x2,…,xm $x_1,x_2,\dots ,x_m$ 的权值叫做回归系数，表达为 θ1,θ2,…,θm ${\theta }_1,{\theta }_2,\dots ,{\theta }_m$。他们的加权和就是你的总得分了。

2.4logistic回归和多元线性回归的区别

logistic回归就是一个线性分类模型，它与线性回归的不同点在于：
表达式不同：
logistic为： p(y=1|x;θ)=δ(θTx)=11+e−θTx=eθTx1+eθTx $p(y=1|x;\theta )=\delta ({\theta }^{T}x)=\frac{1}{1+e^{-{\theta }^{T}x}}=\frac{e^{{\theta }^{T}x}}{1+e^{{\theta }^{T}x}}$, p(y=0|x;θ)=1−p(y=1|x;θ)=11+eθTx $p(y=0|x;\theta )=1-p(y=1|x;\theta )=\frac{1}{1+e^{{\theta }^{T}x}}$
而线性回归方程为： f(x)=θTx $f(x)={\theta }^{T}x$
注： θ是一个权重矩阵 $\theta 是一个权重矩阵$
所以说，Logistic Regression 就是一个被logistic方程归一化后的线性回归，仅此而已。

2.5代价函数的确立

LogisticRegression最基本的学习算法是最大似然。什么是最大似然，可以看看这篇博文“ 从最大似然到EM算法浅解”。
假设我们有n个独立的训练样本 。那每一个观察到的样本 (xi,yi) $(x_i,y_i)$出现的概率是：

p(yi,xi)=p(yi=1|xi)yi(1−p(yi=1|xi)1−yi) $$p(y_i,x_i)=p(y_i=1|x_i{)}^{y_i}(1-p(y_i=1|x_i{)}^{1-y_i})$$

上面为什么是这样呢？依据最大似然函数的定义：已知某个参数能使这个样本出现的概率最大。当 y=1 $y=1$的时候，后面那一项为0，只剩下x属于1类的概率，当 y=0 $y=0$的时候，前面那一项为0，只剩下后面那个 x∈0 $x\in 0$的概率（1减去x属于1的概率）。
所以不管y是0还是1，上面得到的数，都是(x, y)出现的概率。
上述表述的知识针对单个样本，那么， 针对整个样本集，也就是n个独立的样本出现的似然函数（因为样本与样本之间彼此独立，所以n个样本出现的概率就是所有样本概率的乘积 ∏ni=1p(xi) $\underset{i=1}{\overset{n}{\prod }}p(x_i)$为：

L(θ)=∏i=1np(yi=1|xi)yi(1−p(yi=1|xi)1−yi) $$L(\theta )=\underset{i=1}{\overset{n}{\prod }}p(y_i=1|x_i{)}^{y_i}(1-p(y_i=1|x_i{)}^{1-y_i})$$

最大似然法就是求解模型中使得似然函数最大的系数θ值。这个 L(θ) $L(\theta )$ 就是我们的代价函数（cost function）了。
我们先变换下L(θ)：取自然对数，然后化简,最后对 θ $\theta$ 求导
为了简化计算，首先对 θ $\theta$ 进行化简
logL(θ)=log∏ni=1p(yi=1|xi)yi(1−p(yi=1|xi)1−yi) $\log L(\theta )=\log \underset{i=1}{\overset{n}{\prod }}p(y_i=1|x_i{)}^{y_i}(1-p(y_i=1|x_i{)}^{1-y_i})$
=∑ni=1yilogp(yi=1|xi)+(1−yi)log(1−p(yi=1|xi)) $=\underset{i=1}{\overset{n}{\sum }}{y}_i\log p(y_i=1|x_i)+(1-y_i)\log (1-p(y_i=1|x_i))$
=∑ni=1yilogp(yi=1|xi)1−p(yi=1|xi)+∑ni=1log(1−p(yi=1|xi)) $=\underset{i=1}{\overset{n}{\sum }}{y}_i\log \frac{p(y_i=1|x_i)}{1-p(y_i=1|x_i)}+\underset{i=1}{\overset{n}{\sum }}\log (1-p(y_i=1|x_i))$
=∑ni=1yi(θ0+θ1x1+⋯+θmxm)+∑ni=1log(1−p(yi=1|xi)) $=\underset{i=1}{\overset{n}{\sum }}{y}_i({\theta }_0+{\theta }_1{x}_1+\cdots +{\theta }_m{x}_m)+\underset{i=1}{\overset{n}{\sum }}\log (1-p(y_i=1|x_i))$
=∑ni=1yi(θTxi)−∑ni=1log(1+eθTxi) $=\underset{i=1}{\overset{n}{\sum }}{y}_i({\theta }^T{x}_i)-\underset{i=1}{\overset{n}{\sum }}\log (1+e^{{\theta }^T{x}_i})$
上述式子，对 θ $\theta$ 求导得：
αlogL(θ)αθ=∑ni=1yixi−∑ni=1eθTx+i1+eθTxi=∑ni=1(yi−δ(θTxi))xi $\frac{\alpha \log L(\theta )}{\alpha \theta }=\underset{i=1}{\overset{n}{\sum }}{y}_i{x}_i-\underset{i=1}{\overset{n}{\sum }}\frac{e^{{\theta }^{T}x+i}}{1+e^{{\theta }^T{x}_i}}=\underset{i=1}{\overset{n}{\sum }}(y_i-\delta ({\theta }^T{x}_i))x_i$

因为无法确定数据的分布，因而，在生活中不太好确定 θ $\theta$

 在线性回归中，在目标值t服从高斯分布的假设下，最大似然解等价于最小化平方误差函数,此时有解析解（误差函数对w求偏导等于0）。
 在逻辑回归二分类问题中，在目标值C服从伯努利分布的假设下，最大似然解等价于最小化交叉熵误差函数，此时没有解析解。
why?

所以，利用迭代法来计算(迭代法生活中太常见了，什么带约束的优化(KKT)，不带约束的优化，凸函数啊，之类的)，也就是经典的梯度下降(本质就是求偏导)方法。

2.6 优化求解

2.6.1 梯度下降(gradient descent)

Gradient descent 又叫 steepest descent，梯度下降法是一个最优化算法，通常也称为最速下降法。最速下降法是求解无约束优化问题最简单和最古老的方法之一，虽然现已不具有实用性，但是许多有效算法都是以它为基础进行改进和修正而得到的。
最速下降法是用负梯度方向为搜索方向的，最速下降法越接近目标值，步长越小，前进越慢，，在实际过程中，为了使得目标函数能尽快得到收敛，一般是先大步走，然后小步走。梯度下降法的计算过程就是沿梯度下降(负梯度)的方向求解极小值（也可以沿梯度上升方向求解极大值）。
举个例子(参考别人博客的，这个例子很好)狂戳链接：
　首先来看看梯度下降的一个直观的解释。倘若我们现在在一座大山上的某处位置，由于我们不知道怎么下山，于是决定走一步算一步，也就是在每走到一个位置的时候，求解当前位置的梯度，沿着梯度的负方向，也就是沿着当前最陡峭的位置向下走一步，然后继续求解当前位置梯度，即在当前位置沿着最陡峭离地距离最近的位置再往下走一步。这样一步步的走下去，一直走到我们觉得已经到了山脚为止。很显然，如若我们一直这样走下去，有可能我们不能走到山脚，而是到了某一个局部的山峰低处，比如说可能走到了鞍部。

　从上面的例子可以看出，梯度下降不一定能够找到全局的最优解，有可能是一个局部最优解。当然，如果损失函数是凸函数，梯度下降法得到的解就一定是全局最优解（凸函数简直太神奇了，有木有啊，倘若我们在以后的过程中，能证明我们优化的问题是一个凸问题，那么我们最终求得的参数就是全局最优的，，如果对凸函数有兴趣的，可以关注我的另一篇博文，狂戳( http://blog.csdn.net/u012513618/article/details/78392940）。
要找最小值，只需要每一步都往下走（也就是沿着与梯度相反的方向，梯度的方向就是数据变化最大的方向），然后不断的走，便能走到最小值的地方。

但同时也需要更快的到达最小值，所以需要每一步都找下坡最快的地方，也就是每一步走某个方向，都比走其他方法，要离最小值更近。而这个下坡最快的方向，就是梯度的负方向了。
对logistic Regression来说，梯度下降算法新鲜出炉，如下：
θt+1=θt−α∂L(θ)∂θ ${\theta }^{t+1}={\theta }^t-\alpha \frac{\mathrm{\partial }L(\theta )}{\mathrm{\partial }\theta }$
=θt−α(∑ni=1eθTxi1+eθTxi) $={\theta }^t-\alpha (\underset{i=1}{\overset{n}{\sum }}\frac{e^{{\theta }^T{x}_i}}{1+e^{{\theta }^T{x}_i}})$
=θt−α∑ni=1(yi−δ(θTxi))xi $={\theta }^t-\alpha \underset{i=1}{\overset{n}{\sum }}(y_i-\delta ({\theta }^T{x}_i))x_i$
其中，参数α是学习率。
梯度下降算法的伪代码如下：

初始化回归系数为1
重复下面步骤直到收敛
{ 计算整个数据集的梯度 使用alpha x gradient来更新回归系数 }
返回回归系数值

梯度下降法

def GradDescent(dataMatrix, classLabels,maxcycle):
    m,n = shape(dataMatrix)
    alpha = 0.01
    weights = ones((n,1))##array类型
    for i in range(maxCycle):
        h = sigmid(dataMatrix*weights)
        error = classLabels-h
        weights = weights - alpha *dataMatrix.transpose()*error##梯度下降
    return weights            

2.6.2 随机梯度下降SGD (stochastic gradient descent)

通过上面的代码可以看出：
梯度下降算法在每次更新回归系数的时候都需要遍历整个数据集（计算整个数据集的回归误差），那么倘若遇到有数十亿样本和成千上万的特征时，它的计算复杂度太高（无论是从时间还是空间上），另外，这种方法无法解决有实时需求的问题(如增量学习)。
改进的方法一次仅用一个样本点（的回归误差）来更新回归系数。这个方法叫随机梯度下降算法，它是一种增量学习算法。外循环控制epoch次数，内循环用来更新每个样本的回归系数。我们只需要用新样本x来更新已有分类器h的参数即可），增量学习是一种在线学习算法，每次处理的数据集的叫“批处理”(mini-batch)，分多个批次进行计算。
随机梯度下降算法的伪代码如下：

初始化回归系数为1
重复下面步骤直到收敛
{ 对数据集中每个样本 计算该样本的梯度 使用alpha x gradient来更新回归系数 }
返回回归系数值

随机梯度下降法代码

def stocGradDscent(dataMatrix, classLabels):
    m,n = shape(dataMatrix)
    alpha = 0.01
    weights = ones(n)##array类型
    for i in range(m):
        h = sigmid(dataMatrix[i]*weights)##dataMatrix是行向量
        error = classLabels[i]-h
        weights = weights - alpha * error*dataMatrix[i]##梯度下降
    return weights            

2.6.3 改进的随机梯度下降

评价一个优化算法的优劣(时间复杂度度和空间复杂度)主要是看它是否收敛，收敛速度是否快？是否可达到局部最优？

上图展示了随机梯度下降算法在200次迭代中，我们的数据库有100个二维样本，每个样本都对系数调整一次，所以共有200*100=20000次调整）三个回归系数的变化过程。其中系数X2经过50次迭代就达到了稳定值。但系数 X1 $X_1$ 和 X0 $X_0$ 到100次迭代后稳定。而且可恨的是系数 X1 $X_1$ 和 X2 $X_2$ 还在很调皮的周期波动，迭代次数很大了，心还停不下来。
产生这个现象的原因是存在一些无法正确分类的样本点，也就是我们的数据集 线性不可分，但logistic regression是线性分类模型，对非线性可分情况无能为力。
然而我们的优化程序并没能意识到这些不正常的样本点，仍然一视同仁的对待，并没有调整系数去减少对这些样本的分类误差，从而导致了在每次迭代时引发系数的剧烈改变。

那么如何使得算法能避免来回波动，从而快速稳定和收敛到某个值呢？

对随机梯度下降算法，提出两处改进来避免上述的波动：
1）在每次迭代时，调整更新步长alpha的值。
随着迭代的进行，alpha越来越小，这会缓解系数的高频波动（也就是每次迭代系数改变得太大，跳的跨度太大）。当然了，为了避免alpha随着迭代不断减小到接近于0（这时候，系数几乎没有调整，那么迭代也没有意义了），我们约束alpha一定大于一个稍微大点的常数项，具体见代码。
2）每次迭代，改变样本的优化顺序。
也就是随机选择样本来更新回归系数。这样做可以减少周期性的波动，因为样本顺序的改变，使得每次迭代不再形成周期性。

改进的随机梯度下降算法的伪代码如下：

初始化回归系数为1
重复下面步骤直到收敛{ 对随机遍历的数据集中的每个样本 随着迭代的逐渐进行，减小alpha的值 计算该样本的梯度 使用alpha x gradient来更新回归系数 }
返回回归系数值

改进的随机梯度下降算法代码

def ImpstocGradDscent(dataMatrix, classLabels, numIter =150):
    m,n=shape(dataMatrix)
    weights = ones(n)
    for j in range(numIter):
        for i in range(m):
            alpha = 4/(1.0+j+i)+0.01
            randIndex = int(random.uniform(0,len(dataIndex)))
            h = sigmoid(sum(dataMatrix[randIndex]*weights))
            error = classLabels[randIndex] -h
            weights = weights - alpha *error *dataMatrix[randIndex]
            del(dataIndex[randIndex])
   return weights 

比较原始的随机梯度下降和改进后的梯度下降，可以看到两点不同：
1） 系数不再出现周期性波动。
2）系数可以很快的稳定下来，也就是快速收敛。
这里只迭代了20次就收敛了。而上面的随机梯度下降需要迭代200次才能稳定。

3 逻辑回归的适用条件

(1) 因变量为二分类的分类变量或某事件的发生率，并且是数值型变量。但是需要注意，重复计数现象指标不适用于Logistic回归。
(2)残差和因变量都要服从二项分布。二项分布对应的是分类变量，所以不是正态分布，也不是用最小二乘法，而是用最大似然法来解决方程估计和检验问题。
(3) 自变量和Logistic概率是线性关系
(4) 各观测对象间相互独立。
(5)逻辑回归不适用于不均衡样本数据
针对样本不均衡时，有如下处理方式

情形1：正样本>>负样本 且正负样本数量级都很大，采用down-sampling
情形2：正样本>>负样本 量不大
      1).采集更多的数据-> 变成情形1
      2）.over-sampling(图像数据中的旋转和镜像)
      3)修改loss function【有风险，可能会破坏数据本身的凹凸性】    

注：
关于凹凸性，可以见我另外一篇文章http://blog.csdn.net/u012513618/article/details/78392940

代码源码：

# coding=utf-8
# __author__=Eshter Yuu
#无需言，做自己


''' sigmoid函数和logistic回归分类器 梯度下降算法 数据中的缺失项处理 '''


''' 优点：计算代价不高，易于理解和实现 缺点：容易欠拟合，分类精度可能不高 适合数据类型，数值型和标称型数据 '''

''' 梯度上升的伪代码 每个回归系数初始化为1 重复R次： 计算整个数据集的梯度 使用alpha*gradient更新回归系数的向量 返回回归系数 '''
import numpy as np
import matplotlib.pyplot as plt

##logistic回归梯度上升优化算法
def loadDataSet():
    dataMat =[]; labelMat =[]
    fr =open('testSet.txt')
    for line in fr.readlines():
        linArr = line.strip().split()
        dataMat.append([1.0,float(linArr[0]), float(linArr[1])])
        labelMat.append(int(linArr[2]))
    return dataMat, labelMat

def sigmoid(inX):
    return 1.0 / (1+np.exp(-inX))


def gradGradient(dataMatIn, classLabels):
    dataMatrix = np.mat(dataMatIn)
    labelMat = np.mat(classLabels).transpose()#转化为矩阵
    m,n = np.shape(dataMatrix)
    alpha =0.001
    maxCycles = 500
    weights = np.ones((n,1))
    for k in range(maxCycles):
        h = sigmoid(dataMatrix * weights)
        error = (labelMat - h)
        weights = weights +alpha * dataMatrix.transpose()*error
    return weights


###画出数据集和logistic回归最佳拟合直线的函数
def plotBestFit(weights):
    dataMat, labelMat = loadDataSet()
    dataArr =np.array(dataMat)
    n = np.shape(dataArr)[0]
    xcord1 = [];ycord1=[]
    xcord2 =[];ycord2 =[]
    for i in range(n):
        if int(labelMat[i]) ==1:
            xcord1.append(dataArr[i,1]); ycord1.append(dataArr[i,2])
        else:
            xcord2.append(dataArr[i,1]); ycord2.append(dataArr[i,2])
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(xcord1,ycord1,s =30, c ='red', marker ='s')
    ax.scatter(xcord2,ycord2, s =30, c ='green')
    x = np.arange(-3.0, 3.0, 0.1)
    y = (-weights[0]-weights[1]*x)/weights[2]
    ax.plot(x,y)
    plt.xlabel('X1');plt.ylabel('X2')
    plt.show()



''' 随机梯度上升算法 所有回归系数初始化为1 对数据中每个样本 计算该样本的梯度 使用alpha*gradient更新回归系数值 返回回归系数值 '''
def stoGradAscent0(dataMatrix, classLabels):
    ##array类型a*b就是对应元素相乘
    #np.dot(a,b)就是矩阵的乘法
    m,n = np.shape(dataMatrix)##dataMatrix也是array类型
    alpha =0.01
    weights =np.ones(n)##array类型
    for i in range(m):
        h =sigmoid(sum(dataMatrix[i] * weights))
        error = classLabels[i]-h
        weights = weights + alpha * error * dataMatrix[i]
    return weights
''' 改进的随机梯度上升算法 alpha每次迭代都需要调整 '''
def stoGradAscent1(dataMatrix, classLabels, numIter =150):
    m,n=  np.shape(dataMatrix)
    weights = np.ones(n);
    for j in range(numIter):
        dataIndex = list(range(m))##在python3中range不支持del函数，所以得更换数据类型为list
        for i in range(m):
            alpha = 4/ (1.0+j+i)+0.01
            randIndex = int(np.random.uniform(0, len(dataIndex)))
            h = sigmoid(sum(dataMatrix[randIndex]*weights))
            error = classLabels[randIndex] - h
            weights = weights+ alpha *error* dataMatrix[randIndex]
            del(dataIndex[randIndex])##每次迭代每个样本都迭代到
    return weights

''' 处理数据中的缺失值 1使用可用特征的均值来填补缺失值 2使用特殊值来填补缺失值，如-1 3忽略有缺失的样本值 4使用相似样本的均值添补缺失值 5使用量外机器学习算法预测缺失值 '''


''' logistic回归分类函数 '''
def classfyVector(inX, weights):
    prob = sigmoid(sum(inX *weights))
    if prob > 0.5 : return 1.0
    else: return 0.0

def coliTest():
    frTrain = open('horseColicTraining.txt')
    frTest = open('horseColicTest.txt')
    trainingSet = [];trainingLabel =[]
    for line in frTrain.readlines():
        curLine = line.strip().split('\t')
        lineArr = []
        for i in range(21):
            lineArr.append(float(curLine[i]))
        trainingSet.append(lineArr)
        trainingLabel.append(float(curLine[21]))
    trainWeights = stoGradAscent1(np.array(trainingSet),trainingLabel,500)
    errorCount = 0; numTestVec =0.0
    for line in frTest.readlines():
        numTestVec +=1.0
        curLine = line.strip().split('\t')
        lineArr =[]
        for i in range(21):
            lineArr.append(float(curLine[i]))
        if int(classfyVector(np.array(lineArr), trainWeights)) != int(curLine[21]):
            errorCount +=1
    errorRate = float(errorCount)/ numTestVec
    print(" the error rate of this test is :%f"%errorRate)
    return errorRate


def multiTest():
    numTests =10; errorSum = 0.0
    for k in range(numTests):
        errorSum += coliTest()
    print("after %d iterations the average error rate is:%f:"%(numTests, errorSum/float(numTests)))

#dataArr, labelArr = loadDataSet()
#print(dataArr,'\n',labelArr)
#weights = gradGradient(dataArr,labelArr)
#plotBestFit(weights.getA())


# weights1= stoGradAscent0(np.array(dataArr),labelArr)
# print(weights1)
# plotBestFit(weights1)

# weights2= stoGradAscent1(np.array(dataArr),labelArr)
# print(weights2)
# plotBestFit(weights2)


multiTest()

[参考]：
[1]http://blog.csdn.net/zouxy09/article/details/8537620
[2]http://www.cnblogs.com/pinard/p/5970503.html
[3]https://www.zhihu.com/question/27700702
[4]http://blog.csdn.net/zouxy09/article/details/20319673

Edited by Eshter
Email：[e-mail]
版权归Eshter所有