邏輯回歸

一、邏輯回歸的認識

邏輯回歸是一個用來解決二分類的簡便方法。先來看看邏輯回歸解決二分類的基本思想。

之前寫了線性回歸，現在寫邏輯回歸~都叫回歸，有什么不同呢？

首先，從機器學習的角度說一下。機器學習中，有兩個問題是比較相似的，即預測和分類。通常將模型的輸出是有限的離散值的問題稱為分類問題，而將模型的輸出是連續值的問題稱為預測。不同的兩個問題自然有不同的解決方法，對於預測問題，通常采用回歸分析的方法，比如之前線性回歸對於輸入樣本$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邏輯函數  2 logistic<-function(x){  3   return(1/(1+exp(-x[1])));  4 }  5 
#似然函數  7 ll<-function(theta,x,y){  8   n=ncol(x);  9   m=nrow(x); 10   l<-0; 11   for(i in seq(1,m)){ 12     xi<-matrix(x[i,],ncol=1); 13     hxi<-logistic(t(xi)%*%theta); 14     l<-l+y[i]*log(hxi)+(1-y[i])*log((1-hxi)) 15  } 16   return(-l) #取負值為了使用下面的nlm()函數 17 } 18 
#讀入數據，數據來自Andrew Ng機器學習資料 20 x<-read.table("q1x.txt"); 21 x<-as.matrix(x); 22 x<-cbind(rep(1,99),x); 23 y<-scan("q1y.txt"); 24 y<-matrix(y,ncol=1); 25 n=ncol(x); 26 m=nrow(x); 27 
theta<-rep(0,n); 29 theta<-matrix(theta,ncol=1); 30 result<-nlm(ll,theta,x,y);#nlm()為R運用牛頓法求解函數極小值的函數 31 
est<-result$estimate; 33 

plot(x[1:50,2:3],type="p",col="red",xlim=c(0,8), ylim=c(-4,3)); 36 points(x[51:99,2:3],type="p",col="blue"); 37 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