寫在最前面：Logistic迴歸通過Sigmoid函式接受輸入然後進行預測

首先，介紹一下什麼是Sigmoid函式。大家一定聽過海維賽德階躍函式（Heaviside step function），什麼？沒聽過，好吧，換個名字，單位階躍函式，這個認識吧!

這個函式的問題在於該函式在跳躍點上從0瞬間跳躍到1，這樣處理起來不是很方便

然鵝我們還有另一個具有類似性質的函式Sigmoid函式，計算公式如下：
f(z) = 1/1+e^-z

當x軸的尺度足夠大

import matplotlib.pyplot as plt
import numpy as np


def sigmoid(x):
    # 直接返回sigmoid函式
    return 1. / (1. + np.exp(-x))


def plot_sigmoid():
    # param:起點，終點，間距
    x = np.arange(-30, 30, 0.2)
    y = sigmoid(x)
    plt.plot(x, y)
    plt.show()


if __name__ == '__main__':
    plot_sigmoid()
 

Sigmoid函式的輸入記為z，給出輸入z的公式：

z = w0x0 + w1x1 + w2x2 + w3x3 +…+ wnxn

用向量去理解的話，就是w^T * x，w是一個列向量，大家都知道行向量和列向量相乘得到一個數。其中向量x是分類器的輸入資料，向量w是我們要找的最佳引數。

我們通過梯度上升法尋找最佳迴歸係數。梯度上升法的基本思想是：想要找到某函式的最大值，最好的方法是沿著該函式的梯度放向探尋。
記梯度為▽，則函式f(x,y)的梯度▽f(x,y) =

這個梯度意味著要沿著x的方向移動

然後沿著y的方向移動

這裡引入移動量的大小α，用向量表示梯度上升演算法的迭代公式如下：
w := w + α▽wf(w) 大家肯定會疑惑這裡的 := 什麼意思，其實就是程式語言裡的賦值，程式設計裡一般用 = ，但是這裡用 = 會有歧義。

該公式將會一直迭代，直到達到某個停止條件為止。

好了，基礎介紹到這裡
下面舉個栗子

def loadDataSet(filename):
    dataMat = []; labelMat = []
    fr = open(filename)
    for line in fr.readlines():
        lineArr = line.strip().split()
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
        labelMat.append(int(lineArr[2]))
    return dataMat,labelMat
 

該函式的作用是開啟文字

-0.017612	14.053064	0
-1.395634	4.662541	1
-0.752157	6.538620	0
-1.322371	7.152853	0
0.423363	11.054677	0
0.406704	7.067335	1
0.667394	12.741452	0
-2.460150	6.866805	1
0.569411	9.548755	0
-0.026632	10.427743	0
0.850433	6.920334	1
1.347183	13.175500	0
1.176813	3.167020	1
.
.
.

每行前兩個值分別是x1，x2，第三個值是資料對應的類別標籤。

第二個函式是梯度上升演算法，該函式有兩個引數，第一個引數是dataMatin，是一個二維Numpy陣列，每列分別代表每個不同的特徵，每行則代表每個訓練樣本。dataMatin存放的事100x3的矩陣，alpha是目標移動的步長，maxCycless是迭代次數。
由於txt檔案中包含x1,x2，我們還將x0設為1.0，那麼將三個特徵值作為二維矩陣的列，txt中的100個數據作為行，這就是一個100x3的矩陣，第二個引數是一個1x100的矩陣，列是gradAscent的第二個引數classLabels，為了便於計算，我們使用transpose函式進行轉置，轉成100x1的列向量

labelMat = mat(classLabels).transpose()

然後shape函式計算這個100x3的矩陣的維度，m=100,n=3

m,n = shape(dataMatrix)

步長預設為0.001

alpha = 0.001

迭代次數為500次

maxCycles = 500

ones函式構造了一個全1矩陣，n行1列的矩陣

weights = ones((n,1))

相當於將回歸係數全部初始化為1

迴圈500次，將100x3的矩陣和nx1的矩陣相乘，得到100x1的矩陣，再把這個矩陣作為引數輸入Sigmoid函式，error是標籤減去函式的輸出，這裡計算的是真實類別和預測類別的差值，毫無疑問，h是我們的預測類別。按照該差值的方向調整迴歸係數。

    for k in range(maxCycles):              #heavy on matrix operations
        h = sigmoid(dataMatrix*weights)     #matrix mult
        error = (labelMat - h)              #vector subtraction
        weights = weights + alpha * dataMatrix.transpose()* error #matrix mult

下面是這個函式的完整程式碼

def gradAscent(dataMatIn, classLabels):
    dataMatrix = mat(dataMatIn)             #convert to NumPy matrix
    labelMat = mat(classLabels).transpose() #convert to NumPy matrix
    m,n = shape(dataMatrix)
    alpha = 0.001
    maxCycles = 500
    weights = ones((n,1))
    for k in range(maxCycles):              #heavy on matrix operations
        h = sigmoid(dataMatrix*weights)     #matrix mult
        error = (labelMat - h)              #vector subtraction
        weights = weights + alpha * dataMatrix.transpose()* error #matrix mult
    print(weights)
    return weights

計算出迴歸係數：

[[ 4.12414349]
 [ 0.48007329]
 [-0.6168482 ]]

畫出資料集和Logistic最佳擬合直線的函式

def plotBestFit(weights):
    import matplotlib.pyplot as plt
    dataMat,labelMat=loadDataSet('testSet.txt')
    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(-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()

這個例子用到了100x3的資料集，資料量很小，但是卻需要迭代300次，那麼當樣本達到千萬甚至億的級別的時候，遍歷資料集的成本和代價會非常大。

下面介紹一種改進方法，隨機梯度上升演算法
該演算法每次僅用一個樣本點更新迴歸係數

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

可以看到，隨機梯度上升演算法和梯度上升演算法差別並不大

將每一行和迴歸係數相乘， weights是一個1x3的Numpy陣列，初始化為1

[1. 1. 1.]

然後和每一行的特徵值相乘，得出一個numpy型的float

<class 'numpy.float64'>

最後繪製圖像

今天先到這兒吧，後續看情況更新，大家週末愉快~