機器學習實戰ByMatlab（5）：Logistic Regression

原文出處： Liu_LongPo的專欄（@Liu_LongPo）  

什么叫做回歸呢？舉個例子，我們現在有一些數據點，然后我們打算用一條直線來對這些點進行擬合（該曲線稱為最佳擬合曲線），這個擬合過程就被稱為回歸。

利用Logistic回歸進行分類的主要思想是：

根據現有數據對分類邊界線建立回歸公式，以此進行分類。

這里的”回歸“一詞源于最佳擬合，表示要找到最佳擬合參數集。訓練分類器時的嘴閥就是尋找最佳擬合曲線，使用的是最優化算法。

基于Logistic回歸和Sigmoid函數的分類

優點：計算代價不高，易于理解和實現
缺點：容易欠擬合，分類精度可能不高

使用數據類型：數值型和標稱型數據

Sigmoid函數：

波形如下：

當z為0時，值為0.5，當z增大時，g(z)逼近1，當z減小時，g(z)逼近0

Logistic回歸分類器：

對每一個特征都乘以一個回歸系數，然后把所有結果都相加，再講這個總和代入Sigmoid函數中，從而得到一個范圍在0-1之間的數值。任何大于0.5的數據被分為1，小于0.5的數據被分為0.因此Logistic回歸也被看成是一種概率分布。

分類器的函數形式確定之后，現在的問題就是，如何確定回歸系數？

基于最優化方法的最佳回歸系數確定

Sigmoid函數的輸入記為z，由下面公式得出：

如果采用向量的寫法，則上述公式可以寫成：

其中向量X就是分類器的輸入數據，向量W也就是我們要找到的最佳參數，從而使分類器盡可能更加地精確。接下來將介紹幾種需找最佳參數的方法。

梯度上升法

梯度上升法的基本思想：

要找到某函數的最大值，最好的方法是沿著該函數的梯度方向尋找

這里提一下梯度下降法，這個我們應該會更加熟悉，因為我們在很多代價函數J的優化的時候經常用到它，其基本思想是：

要找到某函數的最小值，最好的方法是沿著該函數的梯度方向的反方向尋找

函數的梯度表示方法如下：

移動方向確定了，移動的大小我們稱之為步長，用α表示，用向量來表示的話，梯度下降算法的迭代公式如下：

該公式已知被迭代執行，直到某個停止條件位置，比如迭代次數達到某個指定值或者算法的誤差小到某個允許的誤差范圍內。

注：梯度下降算法中的迭代公式如下：

Matlab 實現

function weight = gradAscent
%%
clc
close all
clear
%%
data = load('testSet.txt');
[row , col] = size(data);
dataMat = data(:,1:col-1);
dataMat = [ones(row,1) dataMat] ;
labelMat = data(:,col);
alpha = 0.001;
maxCycle = 500;
weight = ones(col,1);
for i = 1:maxCycle
 h = sigmoid((dataMat * weight)');
 error = (labelMat - h');
 weight = weight + alpha * dataMat' * error;
end
figure
scatter(dataMat(find(labelMat(:) == 0),2),dataMat(find(labelMat(:) == 0),3),3);
hold on
scatter(dataMat(find(labelMat(:) == 1),2),dataMat(find(labelMat(:) == 1),3),5);
hold on
x = -3:0.1:3;
y = (-weight(1)-weight(2)*x)/weight(3);
plot(x,y)
hold off
end
function returnVals = sigmoid(inX)
 % 注意這里的sigmoid函數要用點除
 returnVals = 1.0./(1.0+exp(-inX));
end

效圖如下：

由上圖可以看到，回歸效果還是挺不錯的，只有2-4個點分類錯誤。

其實這是的梯度上升算法是批量梯度上升算法，每一次更新參數的時候都要講所有的數據集都代入訓練，效果并不好，下面我們將介紹改進版本：隨機梯度上升算法

隨機梯度上升

梯度上升算法在每次更新回歸系數時都要遍歷整個數據集，該方法在處理100個左右的數據集時尚可，但如果有數十億樣本和成千上萬的特征，那么該方法的復雜度就太高了。一種改進方法是一次僅用一個樣本點來更新回歸系數，該方法就稱為隨機梯度上升法。由于可以在新樣本到來之前對分類器進行增量式更新，因此隨機梯度算法是一個在線學習算法。與”在線學習“相對應，一次處理所有數據被稱作是”批處理“

隨機梯度上升算法可以寫成如下的偽代碼：

所有回歸系數初始化為1
對數據集中的每個樣本
計算該樣本的梯度
使用alpha x gradient 更新回歸系數值
返回回歸系數值

Matlab 代碼實現

function stocGradAscent
%%
%
% Description : LogisticRegression using stocGradAsscent
% Author : Liulongpo
% Time：2015-4-18 10:57:25
%
%%
clc
clear 
close all
%%
data = load('testSet.txt');
[row , col] = size(data);
dataMat = [ones(row,1) data(:,1:col-1)];
alpha = 0.01;
labelMat = data(:,col);
weight = ones(col,1);
for i = 1:row
 h = sigmoid(dataMat(i,:)*weight);
 error = labelMat(i) - h;
 dataMat(i,:)
 weight
 weight = weight + alpha * error * dataMat(i,:)'
end
figure
scatter(dataMat(find(labelMat(:)==0),2),dataMat(find(labelMat(:)==0),3),5);
hold on
scatter(dataMat(find(labelMat(:) == 1),2),dataMat(find(labelMat(:) == 1),3),5);
hold on
x = -3:0.1:3;
y = -(weight(1)+weight(2)*x)/weight(3);
plot(x,y)
hold off

end
function returnVals = sigmoid(inX)
 % 注意這里的sigmoid函數要用點除
 returnVals = 1.0./(1.0+exp(-inX));
end

效果如下：

由上圖可以看出，隨機梯度上升算法分類效果并沒有上面的的梯度上升算法分類效果好。

但是直接比較梯度上升算法和隨機梯度上升算法是不公平的，前者是在整個數據集上迭代500次得到的結果，后者只是迭代了100次。一個判斷算法優劣的可靠方法是看它是否收斂，也就是說求解的參數是否達到了穩定值，是否還會不斷變化。

我們讓隨機梯度上升算法在整個數據集上運行200次，迭代過程中3個參數的變化如下圖：

由上圖可以看到，weight1 最先達到穩定，而weight0和weight2則還需要更多的迭代次數來達到穩定。

此時的分類器跟之前的梯度上升算法的分類效果差不多，如下：

但同時我們也可以看到，三個參數都有不同程度的波動。產生這種現象的原因是存在一些不能被正確分類的樣本點（數據集并非線性可分），在每次迭代的時候都會引起參數的劇烈變化。我們期望算法能避免來回波動，從而收斂到某個值。另外，算法收斂速度也要加快。

改進的隨機梯度上升算法

改進的隨機梯度上升算法的主要兩個改進點如下：

1,每一步調整alpha的值，也就是alpha的值是不嚴格下降的
2.隨機采取樣本來更新回歸參數

matlab代碼如下：

function ImproveStocGradAscent
%%
%
% Description : LogisticRegression using stocGradAsscent
% Author : Liulongpo
% Time：2015-4-18 10:57:25
%
%%
clc
clear 
close all
%%
data = load('testSet.txt');
[row , col] = size(data);
dataMat = [ones(row,1) data(:,1:col-1)];
%alpha = 0.01;
numIter = 20;
labelMat = data(:,col);
weightVal = zeros(3,numIter*row);
weight = ones(col,1);
j = 0;
for k = 1:numIter
 randIndex = randperm(row);
 for i = 1:row
 % 改進點 1
 alpha = 4/(1.0+i+k)+0.01; 
 j = j+1;
 % 改進點 2 
 h = sigmoid(dataMat(randIndex(i),:)*weight);
 % 改進點 2
 error = labelMat(randIndex(i)) - h;
 % 改進點 2
 weight = weight + alpha * error * dataMat(randIndex(i),:)';
 weightVal(1,j) = weight(1);
 weightVal(2,j) = weight(2);
 weightVal(3,j) = weight(3);
 end
end
figure
i = 1:numIter*row;
subplot(3,1,1)
plot(i,weightVal(1,:)),title('weight0')%,axis([0 numIter*row 0.8 7])
j = 1:numIter*row;
subplot(3,1,2)
plot(j,weightVal(2,:)),title('weight1')%,axis([0 numIter*row 0.3 1.2])
k = 1:numIter*row;
subplot(3,1,3)
plot(k,weightVal(3,:)),title('weight2')%,axis([0 numIter*row -1.2 -0.1])
figure
scatter(dataMat(find(labelMat(:)==0),2),dataMat(find(labelMat(:)==0),3),5);
hold on
scatter(dataMat(find(labelMat(:) == 1),2),dataMat(find(labelMat(:) == 1),3),5);
hold on
x = -3:0.1:3;
y = -(weight(1)+weight(2)*x)/weight(3);
plot(x,y,'r')
hold off

end
function returnVals = sigmoid(inX)
 % 注意這里的sigmoid函數要用點除
 returnVals = 1.0./(1.0+exp(-inX));
end

改進點 1 中的alpha會隨著迭代次數的增加不斷減小，但由于代碼中常數0.01的存在，alpha不會減少到0。這樣做是為了保證在多次迭代之后新數據對于參數的更新還有一定的影響。

另一點值得注意的就是，alpha每次減少 1/(k+i) ，k 是迭代次數，i是樣本的下標。所以 alpha 不是嚴格下降的。避免參數的嚴格下降也常見于模擬退火算法等其他優化算法中。

第二個改進的地方如代碼注釋中標記的，這里通過隨機采取樣本來更新回歸參數，這樣能夠減少參數的周期性的波動。

由于alpha的動態變化，我們可以在開始的時候設置比較大的值，代碼中設置0.01，alpha也就是每一次迭代的步長，步長越大，越能夠加快參數的收斂速度。然后ahpha會不嚴格下降，這樣就避免了過擬合現象的發生。至于什么是過擬合已經alpha的選取問題將在下面描述。

迭代20次后效果如下：

由上圖可知，步長alpha動態變化之后，參數的收斂速度加快了很多，這里只是對所有樣本數據集迭代20次，weight0 和 weight2很早就收斂。證明了該算法的優異性。

學習率alpha的選取

首先我們看一下梯度上升算法的核心代碼，如下：

h = sigmoid(dataMat(i,:) * weight);
error = labelMat(i) – h;
weight = weight + alpha * error * dataMat(i,:)’;

第一行做的就是估計分類，第二行計算當前估計與正確分類之間的差error，第三行根據這個error來更新參數weight。

我們在迭代的時候，要做的目標就是最小化 error ，我們令 J 代表 error，令向量 θ 代表weight，則很顯然，J是θ的函數。這里盜用Standfor 機器學習教程的圖，如下：

上圖中的每個箭頭就是每一次迭代的更新步長，第一幅圖我們看到，在最小化 J(θ) 的時候迭代了很多次，這說明什么？說明我們要走很多步才能到達全局最優點，原因就是我們每一步走的距離太短，走得太慢，也就是我們的alpha設置得太小。但是當我們處于最優點附近的時候，這樣有利我們向最優點靠近。

下圖中的每個箭頭也代表走一步，我們可以看到，迭代的時候，每一步都沒有到達最優點，而是在最優點的附近波動。為什么呢？因為步長太大了嘛，明明就在眼前了，半步或者四分之三步就走到了，你卻只能一跨而過，重新再來。但是學習率大的話，在剛開始迭代的時候有利于我們參數的快速收斂，也有利于我們避開局部最小值。

綜合以上兩種情況，我們就應該在開始的時候選取較大的學習率，然后不斷不嚴格減小學習率，這樣才是最優的選擇。

那么，我們開始的學習率應該怎么選取？Andrew Ng 在課程中建議先試試0.01，太大就0.003，太小就0.03….