簡介

Logistic迴歸的目的是尋找一個非線性函式Sigmoid的最佳擬合引數，一般使用梯度上升演算法。
對於有n個屬性的train資料集(X1,X2,...Xn),我們尋找一組迴歸係數(W0,W1,W2,...,Wn)使得函式：
sigmoid(W0+W1*X1+W2*X2+...+Wn*Xn)
最佳擬合train資料集的labels。

演算法描述

初始化迴歸係數（一般初始化為1）
Repeat(N):
    計算資料集的梯度
    使用alpha*gradient更新迴歸係數
返回迴歸係數

Logistic迴歸梯度上升優化演算法

#定義sigmoid函式
def sigmoid
 
(inX):
    return 1.0/(1+exp(-inX))

#梯度上升優化演算法    
def gradAscent(dataMatIn , classLabels , iterNum):
    dataMatrix = mat(dataMatIn)
    labelMat = mat(classLabels).transpose()
    m,n = shape(dataMatrix)
    alpha = 0.001
    weights = ones((n,1))
    for k in range(iterNum):
        h = sigmoid(dataMatrix*weights)
        error = (labelMat - h)
        weights = weights + alpha*dataMatrix.transpose()*error
    return
 
 weights

特別說明一下，dataMatIn是原dataSet去掉最後一列（即labels），然後在第一列加入一列1之後的矩陣。
這是因為我們的迴歸引數中有一個常數項W0，這樣統一之後就可以直接在矩陣dataMatIn上進行操作。

隨機梯度上升演算法

def stocGradAscent(dataMatrix , classLabels):
    m,n = shape(dataMatrix)
    alpha = 0.01
    weights = ones(n)
    for i in range(m):
        h = sigmoid(sum(dataMatrix[i]*weights))
        error = classLabels[i] - h
        weights = weights + alpha * error * dataMatrix[i]
    return
 
 weights

可以看到，隨機梯度上升演算法與梯度上升演算法在程式碼上很相似,但也有一些區別：
1.梯度上升演算法的變數h和誤差error都是向量；而隨機梯度上升演算法中，它們都是數值
2.隨機梯度上升演算法沒有矩陣轉換的過程，所有變數資料都是numpy陣列。

改進的隨機梯度上升演算法

def stocGradAscent1(dataMatrix, classLabels, numIter = 150):
    m,n  = shape(dataMatrix)
    weights = ones(n)
    for j in range(numIter):
        dataIndex = range(m)
        for i in range(m):
            alpha = 4 / (1.0+j+i) + 0.01#alpha每次迭代時需要調整
            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

隨機梯度上升和改進的隨機梯度上升演算法的dataMatrix傳入的都是array型別的引數。

畫出決策邊界

def plotBestFit(weights):
    dataMat , labelMat = loadDataSet()
    dataArr = array(dataMat)
    n = 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 = arange(0,2.0,0.1)
    #g(X) = W0 + W1*X1 +W2*X2,圖中的X,Y分別表示X1,X2
    y = (-weights[0] - weights[1] * x) / weights[2]
    ax.plot(x,y)
    plt.xlabel('X1'); plt.ylabel('X2');
    plt.show()

結果測試

資料（來源於西瓜書）

0.697,0.46,1
0.774,0.376,1
0.634,0.264,1
0.608,0.318,1
0.556,0.215,1
0.403,0.237,1
0.481,0.149,1
0.437,0.211,1
0.666,0.091,0
0.639,0.161,0
0.657,0.198,0
0.593,0.042,0
0.719,0.103,0

測試

dataMat , labelMat = loadDataSet()
weights1 = stocGradAscent1(array(dataMat),labelMat,20000)
plotBestFit(weights1)
weights2 = gradAscent(dataMat , labelMat ， 20000)
plotBestFit(weights2.getA())

執行結果

梯度上升演算法迭代20000次

改進的隨機梯度上升演算法迭代20000次

注：本文主要內容來源於Peter Harrington所著的《機器學習實戰》