逻辑回归

一、逻辑回归的认识

逻辑回归是一个用来解决二分类的简便方法。先来看看逻辑回归解决二分类的基本思想。

之前写了线性回归，现在写逻辑回归~都叫回归，有什么不同呢？

首先，从机器学习的角度说一下。机器学习中，有两个问题是比较相似的，即预测和分类。通常将模型的输出是有限的离散值的问题称为分类问题，而将模型的输出是连续值的问题称为预测。不同的两个问题自然有不同的解决方法，对于预测问题，通常采用回归分析的方法，比如之前线性回归对于输入样本$x$，模型的输出为$y={{h}_{\theta }}(x)={{\theta }^{T}}x$，这里$y$的取值是连续性的。

那么回归分析一般是用来做预测的，可不可以用来做分类呢？试想一下，我们用线性回归来解决二分类问题，如果其输出$y$的取值拿一个阈值$\tau $卡一下，对于输出$y<\tau $的样本分为一类，对于输出$y\ge \tau $的样本分为另一类，这不就好了么？问题是对于线性回归的输出的取值范围是没有大小边界的，那么这个阈值$\tau $怎么取就没法弄了。

与线性回归不同，逻辑回归输出的取值范围是0到1之间，这样的话，选一个阈值$\tau $似乎比较可行了，为了说明怎么选，先看看逻辑函数（也叫sigmod函数）长啥样吧~逻辑函数的定义是这样的$g(z)=\frac{1}{1+{{e}^{-z}}}$，$z\in {{\mathbb{R}}^{1}}$ ，$g(z)\in (0,1)$。它的图形长这样：

现在可以看出我们选的阈值应当是0.5，因为从图形上看，函数值取0.5时对应的自变量取值为0，而整个曲线是关于(0,0.5)中心对称的，这意味这选这个点对两类样本而言是比较公平的（可以这么理解吧~）。也就是说，逻辑函数可以将输入分为两类，第一类$z>0$，$g(z)>0.5$，而第二类$z<0$，$g(z)<0.5$。

一般情况下，我们的输入$x$是多维的啊，而逻辑函数的输入是个标量，所以用逻辑函数对我们的多维样本分类首先要对样本进行一个变换$z={{\theta }^{T}}x$（注意对样本加入常数项截距，即$z={{\theta }_{1}}({{x}_{1}}=1)+{{\theta }_{2}}{{x}_{2}}+...+{{\theta }_{n-1}}{{x}_{n-1}}+{{\theta }_{n}}{{x}_{n}}$，对原来n-1为输入增加取值为1的一维，只是为了方便写成向量形式），变换后的逻辑函数可以写出：$${{h}_{\theta }}(x)=g({{\theta }^{T}}x)=\frac{1}{1+{{e}^{-{{\theta }^{T}}x}}}$$　

那么现在要对样本分成两类了：${{\theta }^{T}}x>0,{{h}_{\theta }}(x)>0.5$的一类和${{\theta }^{T}}x<0,{{h}_{\theta }}(x)<0.5$的一类，而逻辑回归的任务就是根据已有的带类别的样本找到这个$\theta $。（后面第三部分结合似然函数再来分析这个结果）。

二、逻辑回归求解

逻辑回归的基本模型已经有了${{h}_{\theta }}(x)=g({{\theta }^{T}}x)=\frac{1}{1+{{e}^{-{{\theta }^{T}}x}}}$，那么怎么求解呢？好像没说解什么呢^_^，当然是解参数$\theta $了，我们要找个$\theta $使得我们的模型是最优的嘛，参数估计问题~所以，先写出对数似然函数吧！

$\ell (\theta )=\log \prod\limits_{i=1}^{m}{p({{y}^{(i)}}|{{x}^{(i)}};\theta )}$，其中m是样本数目，i是样本编号。感觉怪怪的，突然冒出一个${{y}^{(i)}}$来，而且${p({y^{(i)}}|{x^i},\theta )}$也不知道啊。。。${{y}^{(i)}}$是啥呢，当然是样本类别了，我们不是要分类吗，每个样本本身当然有个类别了，这里是二分类就取0或者1吧。接着来看${p({y^{(i)}}|{x^i},\theta )}$，没法求吧，因为我们还没对${y^{(i)}}|{x^i},\theta$建模啊，模型都没有，概率当然求不了。${{y}^{(i)}}$的取值只有0和1两个，所以它应当是一个伯努利分布，我们需要确定分别取这两个值的概率。回顾一下上面的${{h}_{\theta }}(x)$，其取值范围为（0，1），这不正好嘛，${{y}^{(i)}}$取值的概率本身就是0和1之间的么！好了，干脆直接让$p(y=1|x,\theta )={{h}_{\theta }}(x)$，$p(y=0|x,\theta )=1-{{h}_{\theta }}(x)$得了，这正好符合概率的定义的。那这样做有什么意义呢?回顾上面的内容，当${{h}_{\theta }}(x)$小于0.5时，我们将样本分类为类别0，否则分类为类别1，这样的话取值范围为（0，1）的${{h}_{\theta }}(x)$是不是可以衡量样本属于类别1的概率呢？结合逻辑回归函数的曲线看，$z={{\theta }^{T}}x$比0大的越多，${{h}_{\theta }}(x)$取值越是偏离阈值0.5而离1越近，这意味这这个样本的分类越不模糊，很明确的属于其中一个类，相反${{h}_{\theta }}(x)$取值越是偏离阈值0.5而离0越近，样本越是很明确的属于另一个类别。这样的话，样本归为类别1的概率就是${{h}_{\theta }}(x)$，归类为类别0的概率就是$1-{{h}_{\theta }}(x)$。好了，${y^{(i)}}|{x^i},\theta$的模型好了，即:

$p({{y}^{(i)}}=1|{{x}^{(i)}},\theta )={{h}_{\theta }}({{x}^{(i)}})$，$p({{y}^{(i)}}=0|{{x}^{(i)}},\theta )=1-{{h}_{\theta }}({{x}^{(i)}})$

将他俩写出一个式子

$p({{y}^{(i)}}|{{x}^{(i)}},\theta )={{({{h}_{\theta }}({{x}^{(i)}}))}^{{{y}^{(i)}}}}{{(1-{{h}_{\theta }}({{x}^{(i)}}))}^{1-{{y}^{(i)}}}}$ 

现在，可以接着写似然函数了

$\begin{align}\ell (\theta )&=\log \prod\limits_{i=1}^{m}{p({{y}^{(i)}}|{{x}^{(i)}};\theta )} \\ & =\sum\limits_{i=1}^{m}{\log ({{({{h}_{\theta}}({{x}^{(i)}}))}^{{{y}^{(i)}}}}{{(1-{{h}_{\theta }}({{x}^{(i)}}))}^{1-{{y}^{(i)}}}})} \\& =\sum\limits_{i=1}^{m}{\left( {{y}^{(i)}}\log({{h}_{\theta }}({{x}^{(i)}}))+(1-{{y}^{(i)}})\log((1-{{h}_{\theta}}({{x}^{(i)}})) )\right)} \\\end{align}$         

之前提到有人认为可以简单理解为逻辑回归是对线性回归的结果做了一个逻辑函数的转换，然后用来做二分类，现在有了似然函数就可以理解这种说法不准确了，因为两者求解的目标是不一样的，在逻辑回归的似然函数中，不存在线性回归中优化最小二乘的目标。

现在优化这个似然函数使其取最大值就可以了。

对$\theta $的第$l$个分量${{\theta }_{l}}$求偏导(注意${{h}_{\theta }}({{x}^{(i)}})=g({{\theta }^{T}}{{x}^{(i)}})$且利用逻辑函数求导结果$g{{(z)}^{'}}=g(z)(1-g(z))$)

$\begin{align}\frac{\partial \ell (\theta )}{\partial {{\theta }_{l}}}&=\frac{\partial \sum\limits_{i=1}^{m}{\left( {{y}^{(i)}}\log ({{h}_{\theta }}({{x}^{(i)}}))+(1-{{y}^{(i)}})\log((1-{{h}_{\theta }}({{x}^{(i)}}))) \right)}}{\partial {{\theta }_{l}}} \\& =\sum\limits_{i=1}^{m}{\left( {{y}^{(i)}}\frac{1}{{{h}_{\theta }}({{x}^{(i)}})}+(1-{{y}^{(i)}})\frac{1}{1-{{h}_{\theta }}({{x}^{(i)}})} \right)}\frac{\partial {{h}_{\theta }}({{x}^{(i)}})}{\partial {{\theta }_{l}}} \\& =\sum\limits_{i=1}^{m}{\left( {{y}^{(i)}}\frac{1}{{{h}_{\theta }}({{x}^{(i)}})}-(1-{{y}^{(i)}})\frac{1}{1-{{h}_{\theta }}({{x}^{(i)}})} \right)}\frac{\partial g({{\theta }^{T}}{{x}^{(i)}})}{\partial {{\theta }_{l}}} \\& =\sum\limits_{i=1}^{m}{\left( {{y}^{(i)}}\frac{1}{{{h}_{\theta }}({{x}^{(i)}})}-(1-{{y}^{(i)}})\frac{1}{1-{{h}_{\theta }}({{x}^{(i)}})} \right)}g({{\theta }^{T}}{{x}^{(i)}})(1-g({{\theta }^{T}}{{x}^{(i)}}))\frac{\partial {{\theta }^{T}}{{x}^{(i)}}}{\partial {{\theta }_{l}}} \\ & =\sum\limits_{i=1}^{m}{\left( {{y}^{(i)}}(1-g({{\theta }^{T}}{{x}^{(i)}})-(1-{{y}^{(i)}})g({{\theta }^{T}}{{x}^{(i)}}) \right)}x_{l}^{(i)} \\& =\sum\limits_{i=1}^{m}{({{y}^{(i)}}-}g({{\theta }^{T}}{{x}^{(i)}}))x_{l}^{(i)} \\\end{align}$

要得到解析解，可以令偏导为0，解一下${{\theta }_{l}}$，不过解一下就会发现上面的式子不好解。所以为了求解参数$\theta $，用最优化的方法吧~牛顿法，梯度法之类的。

另外一个问题，这个似然函数能不能求最大值呢？也就是说它万一是个凸函数，就可能有最小值而没最大值了，所以我们需要证明个个函数是个凹函数，这样利用优化方法找到一个局部极大值就是全局最大值了，好了，来证其Hessien矩阵半负定吧~

利用上面对${{\theta }_{l}}$的求导结果，Hessian矩阵第$k$行第$l$列的元素

$\begin{align}{{H}_{kl}}&=\frac{{{\partial }^{2}}\ell (\theta )}{\partial k\partial l} \\& =\frac{\partial \sum\limits_{i=1}^{m}{({{y}^{(i)}}-}g({{\theta }^{T}}{{x}^{(i)}}))x_{l}^{(i)}}{\partial {{\theta }_{k}}} \\ & =\sum\limits_{i=1}^{m}{-g({{\theta }^{T}}{{x}^{(i)}})(1-g({{\theta }^{T}}{{x}^{(i)}}))}x_{k}^{(i)}x_{l}^{(i)} \\\end{align}$

那么$H=\sum\limits_{i=1}^{m}{-g({{\theta }^{T}}{{x}^{(i)}})(1-g({{\theta }^{T}}{{x}^{(i)}}))}{{x}^{(i)}}{{x}^{(i)T}}$

对于任意非零向量$p$有

$\begin{align}{{p}^{T}}Hp&=\sum\limits_{i=1}^{m}{-g({{\theta }^{T}}{{x}^{(i)}})(1-g({{\theta }^{T}}{{x}^{(i)}}))}{{p}^{T}}{{x}^{(i)}}{{x}^{(i)T}}p \\& =-\sum\limits_{i=1}^{m}{g({{\theta }^{T}}{{x}^{(i)}})(1-g({{\theta }^{T}}{{x}^{(i)}}))}{{({{p}^{T}}{{x}^{(i)}})}^{2}} \\\end{align}$

由于$0<g({{\theta }^{T}}{{x}^{(i)}})<1$，所以${{p}^{T}}Hp\le 0$，所以似然函数的Hessien矩阵半负定，它是一个凹函数，这样，利用最优化方法求一个局部极大值就可以了。

三、逻辑回归干了干啥

上面推导了逻辑回归的求解的过程，乱乱的，简单总结一下逻辑回归干了个啥吧~

回顾似然函数$\ell (\theta ) = \log \prod\limits_{i = 1}^m {p({y^{(i)}}|{x^i},\theta )}$，我们的目标要最大化这个东西，也就是要最大化连乘符号里面的每一项$p({{y}^{(i)}}|{{x}^{(i)}},\theta )={{({{h}_{\theta }}({{x}^{(i)}}))}^{{{y}^{(i)}}}}{{(1-{{h}_{\theta }}({{x}^{(i)}}))}^{1-{{y}^{(i)}}}}$，它怎么才能大呢？考虑单个样本$x$如果它对应类别$y=1$，那么$p(y|x,\theta )={{h}_{\theta }}(x)$，所以${{h}_{\theta }}(x)$要比较大才好，而${{h}_{\theta }}(x)=g({{\theta }^{T}}x)=\frac{1}{1+{{e}^{-{{\theta }^{T}}x}}}$，所以${{\theta }^{T}}x$要大于0比较好（结合逻辑函数曲线看看）；相反如果$y=0$，最大化似然函数则要求${{\theta }^{T}}x$尽可能小于零。

所以最大化似然函数的解就是找到一个$\theta $，是得对于类别为1的样本，尽可能满足${{\theta }^{T}}x>0$，而对于类别为0的样本，尽可能满足${{\theta }^{T}}x<0$。换句话说，我们找到的超平面${{\theta }^{T}}x=0$用来对样本分类。可以看出，现在从似然函数分析的结果和我们之前第一部分末尾提出的逻辑回归的目标是一致的。

四、简单实现

#l逻辑函数
logistic<-function(x){
  return(1/(1+exp(-x[1])));
}

#似然函数
ll<-function(theta,x,y){
  n=ncol(x);
  m=nrow(x);
  l<-0;
  for(i in seq(1,m)){
    xi<-matrix(x[i,],ncol=1);
    hxi<-logistic(t(xi)%*%theta);
    l<-l+y[i]*log(hxi)+(1-y[i])*log((1-hxi))
  }
  return(-l) #取负值为了使用下面的nlm()函数
}

#读入数据，数据来自Andrew Ng机器学习资料
x<-read.table("q1x.txt");
x<-as.matrix(x);
x<-cbind(rep(1,99),x);
y<-scan("q1y.txt");
y<-matrix(y,ncol=1);
n=ncol(x);
m=nrow(x);

theta<-rep(0,n);
theta<-matrix(theta,ncol=1);
result<-nlm(ll,theta,x,y);#nlm()为R运用牛顿法求解函数极小值的函数

est<-result$estimate;


plot(x[1:50,2:3],type="p",col="red",xlim=c(0,8), ylim=c(-4,3));
points(x[51:99,2:3],type="p",col="blue");
abline(-est[1]/est[3],-est[2]/est[3]);

 参考资料：

[1]Andrew Ng机器学习视频、讲义：http://cs229.stanford.edu/

[2]《R语言实战》 /(美)科巴科弗(Kabacoff,R.I.)著；高涛，肖楠，陈钢译.人民邮电出版社，2013.1

[3]回归概念学习：http://blog.csdn.net/viewcode/article/details/8794401