05機器學習實戰之Logistic 迴歸

Logistic 迴歸 概述

Logistic 迴歸 或者叫邏輯迴歸 雖然名字有迴歸，可是它是用來作分類的。其主要思想是: 根據現有數據對分類邊界線(Decision Boundary)創建迴歸公式，以此進行分類。

須知概念

Sigmoid 函數

迴歸 概念

假設如今有一些數據點，咱們用一條直線對這些點進行擬合（這條直線稱爲最佳擬合直線），這個擬合的過程就叫作迴歸。進而能夠獲得對這些點的擬合直線方程，那麼咱們根據這個迴歸方程，怎麼進行分類呢？請看下面。

二值型輸出分類函數

咱們想要的函數應該是: 能接受全部的輸入而後預測出類別。例如，在兩個類的狀況下，上述函數輸出 0 或 1.或許你以前接觸過具備這種性質的函數，該函數稱爲 海維塞得階躍函數(Heaviside step function)，或者直接稱爲 單位階躍函數。然而，海維塞得階躍函數的問題在於: 該函數在跳躍點上從 0 瞬間跳躍到 1，這個瞬間跳躍過程有時很難處理。幸虧，另外一個函數也有相似的性質（能夠輸出 0 或者 1 的性質），且數學上更易處理，這就是 Sigmoid 函數。 Sigmoid 函數具體的計算公式以下:

下圖給出了 Sigmoid 函數在不一樣座標尺度下的兩條曲線圖。當 x 爲 0 時，Sigmoid 函數值爲 0.5 。隨着 x 的增大，對應的 Sigmoid 值將逼近於 1 ; 而隨着 x 的減少， Sigmoid 值將逼近於 0 。若是橫座標刻度足夠大， Sigmoid 函數看起來很像一個階躍函數。

所以，爲了實現 Logistic 迴歸分類器，咱們能夠在每一個特徵上都乘以一個迴歸係數（以下公式所示），而後把全部結果值相加，將這個總和代入 Sigmoid 函數中，進而獲得一個範圍在 0~1 之間的數值。任何大於 0.5 的數據被分入 1 類，小於 0.5 即被納入 0 類。因此，Logistic 迴歸也是一種機率估計，好比這裏Sigmoid 函數得出的值爲0.5，能夠理解爲給定數據和參數，數據被分入 1 類的機率爲0.5。想對Sigmoid 函數有更多瞭解，能夠點開此連接跟此函數互動。

基於最優化方法的迴歸係數肯定

Sigmoid 函數的輸入記爲 z ，由下面公式獲得:

若是採用向量的寫法，上述公式能夠寫成  ，它表示將這兩個數值向量對應元素相乘而後所有加起來即獲得 z 值。其中的向量 x 是分類器的輸入數據，向量 w 也就是咱們要找到的最佳參數（係數），從而使得分類器儘量地精確。爲了尋找該最佳參數，須要用到最優化理論的一些知識。咱們這裏使用的是——梯度上升法（Gradient Ascent）。

梯度上升法

梯度上升法的思想

要找到某函數的最大值，最好的方法是沿着該函數的梯度方向探尋。若是梯度記爲 ▽ ，則函數 f(x, y) 的梯度由下式表示:

這個梯度意味着要沿 x 的方向移動  ，沿 y 的方向移動  。其中，函數f(x, y) 必需要在待計算的點上有定義而且可微。下圖是一個具體的例子。

上圖展現的，梯度上升算法到達每一個點後都會從新估計移動的方向。從 P0 開始，計算完該點的梯度，函數就根據梯度移動到下一點 P1。在 P1 點，梯度再次被從新計算，並沿着新的梯度方向移動到 P2 。如此循環迭代，直到知足中止條件。迭代過程當中，梯度算子老是保證咱們能選取到最佳的移動方向。

上圖中的梯度上升算法沿梯度方向移動了一步。能夠看到，梯度算子老是指向函數值增加最快的方向。這裏所說的是移動方向，而未提到移動量的大小。該量值稱爲步長，記做 α 。用向量來表示的話，梯度上升算法的迭代公式以下:

該公式將一直被迭代執行，直至達到某個中止條件爲止，好比迭代次數達到某個指定值或者算法達到某個能夠容許的偏差範圍。

介紹一下幾個相關的概念：

例如：y = w0 + w1x1 + w2x2 + ... + wnxn
梯度：參考上圖的例子，二維圖像，x方向表明第一個係數，也就是 w1，y方向表明第二個係數也就是 w2，這樣的向量就是梯度。
α：上面的梯度算法的迭代公式中的阿爾法，這個表明的是移動步長（step length）。移動步長會影響最終結果的擬合程度，最好的方法就是隨着迭代次數更改移動步長。
步長通俗的理解，100米，若是我一步走10米，我須要走10步；若是一步走20米，我只須要走5步。這裏的一步走多少米就是步長的意思。
▽f(w)：表明沿着梯度變化的方向。

問：有人會好奇爲何有些書籍上說的是梯度降低法（Gradient Decent）?

答： 其實這個兩個方法在此狀況下本質上是相同的。關鍵在於代價函數（cost function）或者叫目標函數（objective function）。若是目標函數是損失函數，那就是最小化損失函數來求函數的最小值，就用梯度降低。 若是目標函數是似然函數（Likelihood function），就是要最大化似然函數來求函數的最大值，那就用梯度上升。在邏輯迴歸中， 損失函數和似然函數無非就是互爲正負關係。

只須要在迭代公式中的加法變成減法。所以，對應的公式能夠寫成

局部最優現象 （Local Optima）

上圖表示參數 θ 與偏差函數 J(θ) 的關係圖 (這裏的偏差函數是損失函數，因此咱們要最小化損失函數)，紅色的部分是表示 J(θ) 有着比較高的取值，咱們須要的是，可以讓 J(θ) 的值儘可能的低。也就是深藍色的部分。θ0，θ1 表示 θ 向量的兩個維度（此處的θ0，θ1是x0和x1的係數，也對應的是上文w0和w1）。

可能梯度降低的最終點並不是是全局最小點，多是一個局部最小點，如咱們上圖中的右邊的梯度降低曲線，描述的是最終到達一個局部最小點，這是咱們從新選擇了一個初始點獲得的。

看來咱們這個算法將會在很大的程度上被初始點的選擇影響而陷入局部最小點。

Logistic 迴歸 原理

Logistic 迴歸 工做原理

每一個迴歸係數初始化爲 1
重複 R 次:
    計算整個數據集的梯度
    使用 步長 x 梯度 更新迴歸係數的向量
返回迴歸係數

Logistic 迴歸 開發流程

收集數據: 採用任意方法收集數據
準備數據: 因爲須要進行距離計算，所以要求數據類型爲數值型。另外，結構化數據格式則最佳。
分析數據: 採用任意方法對數據進行分析。
訓練算法: 大部分時間將用於訓練，訓練的目的是爲了找到最佳的分類迴歸係數。
測試算法: 一旦訓練步驟完成，分類將會很快。
使用算法: 首先，咱們須要輸入一些數據，並將其轉換成對應的結構化數值；接着，基於訓練好的迴歸係數就能夠對這些數值進行簡單的迴歸計算，斷定它們屬於哪一個類別；在這以後，咱們就能夠在輸出的類別上作一些其餘分析工做。

Logistic 迴歸 算法特色

優勢: 計算代價不高，易於理解和實現。
缺點: 容易欠擬合，分類精度可能不高。
適用數據類型: 數值型和標稱型數據。

附加 方向導數與梯度

Logistic 迴歸 項目案例

項目概述

在一個簡單的數據集上，採用梯度上升法找到 Logistic 迴歸分類器在此數據集上的最佳迴歸係數

開發流程

收集數據: 可使用任何方法
準備數據: 因爲須要進行距離計算，所以要求數據類型爲數值型。另外，結構化數據格式則最佳
分析數據: 畫出決策邊界
訓練算法: 使用梯度上升找到最佳參數
測試算法: 使用 Logistic 迴歸進行分類
使用算法: 對簡單數據集中數據進行分類

收集數據: 可使用任何方法

咱們採用存儲在 TestSet.txt 文本文件中的數據，存儲格式以下:

-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

繪製在圖中，以下圖所示:

準備數據: 因爲須要進行距離計算，所以要求數據類型爲數值型。另外，結構化數據格式則最佳

import numpy as np

data_arr = []
label_arr = []
f = open('D:\\mlInAction\\data\\5.Logistic\\TestSet.txt', 'r')
for line in f.readlines():
    line_arr = line.strip().split()
    # 爲了方便計算，咱們將 X0 的值設爲 1.0 ，也就是在每一行的開頭添加一個 1.0 做爲 X0
    data_arr.append([1.0, np.float(line_arr[0]), np.float(line_arr[1])])
    label_arr.append(int(line_arr[2]))

分析數據: 採用任意方法對數據進行分析，此處不須要

訓練算法: 使用梯度上升找到最佳參數

定義sigmoid階躍函數

def sigmoid(x):
    # 這裏其實很是有必要解釋一下，會出現的錯誤 RuntimeWarning: overflow encountered in exp
    # 這個錯誤在學習階段雖然能夠忽略，可是咱們至少應該知道爲何
    # 這裏是由於咱們輸入的有的 x 實在是過小了，好比 -6000之類的，那麼計算一個數字 np.exp(6000)這個結果太大了，無法表示，因此就溢出了
    # 若是是計算 np.exp（-6000），這樣雖然也會溢出，可是這是下溢，就是表示成零
    # 去網上搜了不少方法，好比 使用bigfloat這個庫（我居然沒有安裝成功，就不嘗試了，反正應該是有用的
    return 1.0 / (1 + np.exp(-x))

Logistic 迴歸梯度上升優化算法

def grad_ascent(data_arr, class_labels):
    """
    梯度上升法，其實就是由於使用了極大似然估計，這個你們有必要去看推導，只看代碼感受不太夠
    :param data_arr: 傳入的就是一個普通的數組，固然你傳入一個二維的ndarray也行
    :param class_labels: class_labels 是類別標籤，它是一個 1*100 的行向量。
                    爲了便於矩陣計算，須要將該行向量轉換爲列向量，作法是將原向量轉置，再將它賦值給label_mat
    :return: 
    """
    # 注意一下，我把原來 data_mat_in 改爲data_arr,由於傳進來的是一個數組，用這個比較不容易搞混
    # turn the data_arr to numpy matrix
    data_mat = np.mat(data_arr)
    # 變成矩陣以後進行轉置
    label_mat = np.mat(class_labels).transpose()
    # m->數據量，樣本數 n->特徵數
    m, n = np.shape(data_mat)
    # 學習率，learning rate
    alpha = 0.001
    # 最大迭代次數，僞裝迭代這麼屢次就能收斂2333
    max_cycles = 500
    # 生成一個長度和特徵數相同的矩陣，此處n爲3 -> [[1],[1],[1]]
    # weights 表明迴歸係數， 此處的 ones((n,1)) 建立一個長度和特徵數相同的矩陣，其中的數所有都是 1
    weights = np.ones((n, 1))
    for k in range(max_cycles):
        # 這裏是點乘  m x 3 dot 3 x 1
        h = sigmoid(data_mat * weights)
        error = label_mat - h
        # 這裏比較建議看一下推導，爲何這麼作能夠，這裏已是求導以後的
        weights = weights + alpha * data_mat.transpose() * error
    return weights

你們看到這兒可能會有一些疑惑，就是，咱們在迭代中更新咱們的迴歸係數，後邊的部分是怎麼計算出來的？爲何會是 alpha * dataMatrix.transpose() * error ?由於這就是咱們所求的梯度，也就是對 f(w) 對 w 求一階導數。具體推導以下:

 

可參考http://blog.csdn.net/achuo/article/details/51160101

畫出數據集和 Logistic 迴歸最佳擬合直線的函數

import matplotlib.pyplot as plt

def plot_best_fit(data_mat, label_mat, weights):
    """
    可視化
    :param weights: 
    :return: 
    """
    data_arr = np.array(data_mat)
    n = np.shape(data_mat)[0]
    x_cord1 = []
    y_cord1 = []
    x_cord2 = []
    y_cord2 = []
    for i in range(n):
        if int(label_mat[i]) == 1:
            x_cord1.append(data_arr[i, 1])
            y_cord1.append(data_arr[i, 2])
        else:
            x_cord2.append(data_arr[i, 1])
            y_cord2.append(data_arr[i, 2])
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(x_cord1, y_cord1, s=30, color='k', marker='^')
    ax.scatter(x_cord2, y_cord2, s=30, color='red', marker='s')
    x = np.arange(-3.0, 3.0, 0.1)
    # print(x)
    y = (-weights[0] - weights[1] * x) / weights[2]
    # type(y)
    y = np.ravel(y)  # y原來是一個二維，須要轉化爲1維
    """
    y的由來，臥槽，是否是沒看懂？
    首先理論上是這個樣子的。
    dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
    w0*x0+w1*x1+w2*x2=f(x)
    x0最開始就設置爲1叻， x2就是咱們畫圖的y值，而f(x)被咱們磨合偏差給算到w0,w1,w2身上去了
    因此： w0+w1*x+w2*y=0 => y = (-w0-w1*x)/w2   
    """
    ax.plot(x, y)
    plt.xlabel('x1')
    plt.ylabel('y1')
    plt.show()

注意

梯度上升算法在每次更新迴歸係數時都須要遍歷整個數據集，該方法在處理 100 個左右的數據集時尚可，但若是有數十億樣本和成千上萬的特徵，那麼該方法的計算複雜度就過高了。一種改進方法是一次僅用一個樣本點來更新迴歸係數，該方法稱爲 隨機梯度上升算法。因爲能夠在新樣本到來時對分類器進行增量式更新，於是隨機梯度上升算法是一個在線學習(online learning)算法。與 「在線學習」 相對應，一次處理全部數據被稱做是 「批處理」 （batch） 。

隨機梯度上升算法能夠寫成以下的僞代碼:

如下是隨機梯度上升算法的實現代碼:全部迴歸係數初始化爲 1
對數據集中每一個樣本
    計算該樣本的梯度
    使用 alpha * gradient 更新迴歸係數值
返回迴歸係數值

def stoc_grad_ascent0(data_mat, class_labels):
    """
    隨機梯度上升，只使用一個樣本點來更新迴歸係數
    :param data_mat: 輸入數據的數據特徵（除去最後一列）,ndarray
    :param class_labels: 輸入數據的類別標籤（最後一列數據）
    :return: 獲得的最佳迴歸係數
    """
    m, n = np.shape(data_mat)
    alpha = 0.01
    weights = np.ones(n)
    for i in range(m):
        # sum(data_mat[i]*weights)爲了求 f(x)的值， f(x)=a1*x1+b2*x2+..+nn*xn,
        # 此處求出的 h 是一個具體的數值，而不是一個矩陣
        h = sigmoid(sum(data_mat[i] * weights))
        error = class_labels[i] - h
        # 仍是和上面同樣，這個先去看推導，再寫程序
        weights = weights + alpha * error * data_mat[i]
    return weights

能夠看到，隨機梯度上升算法與梯度上升算法在代碼上很類似，但也有一些區別: 第一，後者的變量 h 和偏差 error 都是向量，而前者則全是數值；第二，前者沒有矩陣的轉換過程，全部變量的數據類型都是 NumPy 數組。

判斷優化算法優劣的可靠方法是看它是否收斂，也就是說參數是否達到了穩定值，是否還會不斷地變化？下圖展現了隨機梯度上升算法在 200 次迭代過程當中迴歸係數的變化狀況。其中的係數2，也就是 X2 只通過了 50 次迭代就達到了穩定值，但係數 1 和 0 則須要更屢次的迭代。以下圖所示:

針對這個問題，咱們改進了以前的隨機梯度上升算法，以下:

def stoc_grad_ascent1(data_mat, class_labels, num_iter=150):
    """
    改進版的隨機梯度上升，使用隨機的一個樣原本更新迴歸係數
    :param data_mat: 輸入數據的數據特徵（除去最後一列）,ndarray
    :param class_labels: 輸入數據的類別標籤（最後一列數據
    :param num_iter: 迭代次數
    :return: 獲得的最佳迴歸係數
    """
    m, n = np.shape(data_mat)
    weights = np.ones(n)
    for j in range(num_iter):
        # 這裏必需要用list，否則後面的del無法使用
        data_index = list(range(m))
        for i in range(m):
            # i和j的不斷增大，致使alpha的值不斷減小，可是不爲0
            alpha = 4 / (1.0 + j + i) + 0.01
            # 隨機產生一個 0～len()之間的一個值
            # random.uniform(x, y) 方法將隨機生成下一個實數，它在[x,y]範圍內,x是這個範圍內的最小值，y是這個範圍內的最大值。
            rand_index = int(np.random.uniform(0, len(data_index)))
            h = sigmoid(np.sum(data_mat[data_index[rand_index]] * weights))
            error = class_labels[data_index[rand_index]] - h
            weights = weights + alpha * error * data_mat[data_index[rand_index]]
            del(data_index[rand_index])
    return weights

In [54]:

import numpy as np

data_arr = []
label_arr = []
f = open('D:\\mlInAction\\data\\5.Logistic\\TestSet.txt', 'r')
for line in f.readlines():
    line_arr = line.strip().split()
    # 爲了方便計算，咱們將 X0 的值設爲 1.0 ，也就是在每一行的開頭添加一個 1.0 做爲 X0
    data_arr.append([1.0, np.float(line_arr[0]), np.float(line_arr[1])])
    label_arr.append(int(line_arr[2]))

In [55]:

data_arr

Out[55]:

[[1.0, -0.017612, 14.053064],
 [1.0, -1.395634, 4.662541],
 [1.0, -0.752157, 6.53862],
 [1.0, -1.322371, 7.152853],
 [1.0, 0.423363, 11.054677],
 [1.0, 0.406704, 7.067335],
 [1.0, 0.667394, 12.741452],
 [1.0, -2.46015, 6.866805],
 [1.0, 0.569411, 9.548755],
 [1.0, -0.026632, 10.427743],
 [1.0, 0.850433, 6.920334],
 [1.0, 1.347183, 13.1755],
 [1.0, 1.176813, 3.16702],
 [1.0, -1.781871, 9.097953],
 [1.0, -0.566606, 5.749003],
 [1.0, 0.931635, 1.589505],
 [1.0, -0.024205, 6.151823],
 [1.0, -0.036453, 2.690988],
 [1.0, -0.196949, 0.444165],
 [1.0, 1.014459, 5.754399],
 [1.0, 1.985298, 3.230619],
 [1.0, -1.693453, -0.55754],
 [1.0, -0.576525, 11.778922],
 [1.0, -0.346811, -1.67873],
 [1.0, -2.124484, 2.672471],
 [1.0, 1.217916, 9.597015],
 [1.0, -0.733928, 9.098687],
 [1.0, -3.642001, -1.618087],
 [1.0, 0.315985, 3.523953],
 [1.0, 1.416614, 9.619232],
 [1.0, -0.386323, 3.989286],
 [1.0, 0.556921, 8.294984],
 [1.0, 1.224863, 11.58736],
 [1.0, -1.347803, -2.406051],
 [1.0, 1.196604, 4.951851],
 [1.0, 0.275221, 9.543647],
 [1.0, 0.470575, 9.332488],
 [1.0, -1.889567, 9.542662],
 [1.0, -1.527893, 12.150579],
 [1.0, -1.185247, 11.309318],
 [1.0, -0.445678, 3.297303],
 [1.0, 1.042222, 6.105155],
 [1.0, -0.618787, 10.320986],
 [1.0, 1.152083, 0.548467],
 [1.0, 0.828534, 2.676045],
 [1.0, -1.237728, 10.549033],
 [1.0, -0.683565, -2.166125],
 [1.0, 0.229456, 5.921938],
 [1.0, -0.959885, 11.555336],
 [1.0, 0.492911, 10.993324],
 [1.0, 0.184992, 8.721488],
 [1.0, -0.355715, 10.325976],
 [1.0, -0.397822, 8.058397],
 [1.0, 0.824839, 13.730343],
 [1.0, 1.507278, 5.027866],
 [1.0, 0.099671, 6.835839],
 [1.0, -0.344008, 10.717485],
 [1.0, 1.785928, 7.718645],
 [1.0, -0.918801, 11.560217],
 [1.0, -0.364009, 4.7473],
 [1.0, -0.841722, 4.119083],
 [1.0, 0.490426, 1.960539],
 [1.0, -0.007194, 9.075792],
 [1.0, 0.356107, 12.447863],
 [1.0, 0.342578, 12.281162],
 [1.0, -0.810823, -1.466018],
 [1.0, 2.530777, 6.476801],
 [1.0, 1.296683, 11.607559],
 [1.0, 0.475487, 12.040035],
 [1.0, -0.783277, 11.009725],
 [1.0, 0.074798, 11.02365],
 [1.0, -1.337472, 0.468339],
 [1.0, -0.102781, 13.763651],
 [1.0, -0.147324, 2.874846],
 [1.0, 0.518389, 9.887035],
 [1.0, 1.015399, 7.571882],
 [1.0, -1.658086, -0.027255],
 [1.0, 1.319944, 2.171228],
 [1.0, 2.056216, 5.019981],
 [1.0, -0.851633, 4.375691],
 [1.0, -1.510047, 6.061992],
 [1.0, -1.076637, -3.181888],
 [1.0, 1.821096, 10.28399],
 [1.0, 3.01015, 8.401766],
 [1.0, -1.099458, 1.688274],
 [1.0, -0.834872, -1.733869],
 [1.0, -0.846637, 3.849075],
 [1.0, 1.400102, 12.628781],
 [1.0, 1.752842, 5.468166],
 [1.0, 0.078557, 0.059736],
 [1.0, 0.089392, -0.7153],
 [1.0, 1.825662, 12.693808],
 [1.0, 0.197445, 9.744638],
 [1.0, 0.126117, 0.922311],
 [1.0, -0.679797, 1.22053],
 [1.0, 0.677983, 2.556666],
 [1.0, 0.761349, 10.693862],
 [1.0, -2.168791, 0.143632],
 [1.0, 1.38861, 9.341997],
 [1.0, 0.317029, 14.739025]]

In [56]:

label_arr

Out[56]:

[0,
 1,
 0,
 0,
 0,
 1,
 0,
 1,
 0,
 0,
 1,
 0,
 1,
 0,
 1,
 1,
 1,
 1,
 1,
 1,
 1,
 1,
 0,
 1,
 1,
 0,
 0,
 1,
 1,
 0,
 1,
 1,
 0,
 1,
 1,
 0,
 0,
 0,
 0,
 0,
 1,
 1,
 0,
 1,
 1,
 0,
 1,
 1,
 0,
 0,
 0,
 0,
 0,
 0,
 1,
 1,
 0,
 1,
 0,
 1,
 1,
 1,
 0,
 0,
 0,
 1,
 1,
 0,
 0,
 0,
 0,
 1,
 0,
 1,
 0,
 0,
 1,
 1,
 1,
 1,
 0,
 1,
 0,
 1,
 1,
 1,
 1,
 0,
 1,
 1,
 1,
 0,
 0,
 1,
 1,
 1,
 0,
 1,
 0,
 0]

In [57]:

def sigmoid(x):
    # 這裏其實很是有必要解釋一下，會出現的錯誤 RuntimeWarning: overflow encountered in exp
    # 這個錯誤在學習階段雖然能夠忽略，可是咱們至少應該知道爲何
    # 這裏是由於咱們輸入的有的 x 實在是過小了，好比 -6000之類的，那麼計算一個數字 np.exp(6000)這個結果太大了，無法表示，因此就溢出了
    # 若是是計算 np.exp（-6000），這樣雖然也會溢出，可是這是下溢，就是表示成零
    # 去網上搜了不少方法，好比 使用bigfloat這個庫（我居然沒有安裝成功，就不嘗試了，反正應該是有用的
    return 1.0 / (1 + np.exp(-x))

In [58]:

def grad_ascent(data_arr, class_labels):
    """
    梯度上升法，其實就是由於使用了極大似然估計，這個你們有必要去看推導，只看代碼感受不太夠
    :param data_arr: 傳入的就是一個普通的數組，固然你傳入一個二維的ndarray也行
    :param class_labels: class_labels 是類別標籤，它是一個 1*100 的行向量。
                    爲了便於矩陣計算，須要將該行向量轉換爲列向量，作法是將原向量轉置，再將它賦值給label_mat
    :return: 
    """
    # 注意一下，我把原來 data_mat_in 改爲data_arr,由於傳進來的是一個數組，用這個比較不容易搞混
    # turn the data_arr to numpy matrix
    data_mat = np.mat(data_arr)
    # 變成矩陣以後進行轉置
    label_mat = np.mat(class_labels).transpose()
    # m->數據量，樣本數 n->特徵數
    m, n = np.shape(data_mat)
    # 學習率，learning rate
    alpha = 0.001
    # 最大迭代次數，僞裝迭代這麼屢次就能收斂2333
    max_cycles = 500
    # 生成一個長度和特徵數相同的矩陣，此處n爲3 -> [[1],[1],[1]]
    # weights 表明迴歸係數， 此處的 ones((n,1)) 建立一個長度和特徵數相同的矩陣，其中的數所有都是 1
    weights = np.ones((n, 1))
    for k in range(max_cycles):
        # 這裏是點乘  m x 3 dot 3 x 1
        h = sigmoid(data_mat * weights)
        error = label_mat - h
        # 這裏比較建議看一下推導，爲何這麼作能夠，這裏已是求導以後的
        weights = weights + alpha * data_mat.transpose() * error
    return weights

In [69]:

weights = grad_ascent(data_arr, label_arr)
weights

Out[69]:

matrix([[ 4.12414349],
        [ 0.48007329],
        [-0.6168482 ]])

In [76]:

import matplotlib.pyplot as plt


def plot_best_fit(data_mat, label_mat, weights):
    """
    可視化
    :param weights: 
    :return: 
    """
    data_arr = np.array(data_mat)
    n = np.shape(data_mat)[0]
    x_cord1 = []
    y_cord1 = []
    x_cord2 = []
    y_cord2 = []
    for i in range(n):
        if int(label_mat[i]) == 1:
            x_cord1.append(data_arr[i, 1])
            y_cord1.append(data_arr[i, 2])
        else:
            x_cord2.append(data_arr[i, 1])
            y_cord2.append(data_arr[i, 2])
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(x_cord1, y_cord1, s=30, color='k', marker='^')
    ax.scatter(x_cord2, y_cord2, s=30, color='red', marker='s')
    x = np.arange(-3.0, 3.0, 0.1)
    # print(x)
    y = (-weights[0] - weights[1] * x) / weights[2]
    # type(y)
    y = np.ravel(y)  # y原來是一個二維，須要轉化爲1維
    """
    y的由來，臥槽，是否是沒看懂？
    首先理論上是這個樣子的。
    dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
    w0*x0+w1*x1+w2*x2=f(x)
    x0最開始就設置爲1叻， x2就是咱們畫圖的y值，而f(x)被咱們磨合偏差給算到w0,w1,w2身上去了
    因此： w0+w1*x+w2*y=0 => y = (-w0-w1*x)/w2   
    """
    ax.plot(x, y)
    plt.xlabel('x1')
    plt.ylabel('y1')
    plt.show()

In [77]:

plot_best_fit(data_arr, label_arr, weights)

 

In [72]:

def stoc_grad_ascent0(data_mat, class_labels):
    """
    隨機梯度上升，只使用一個樣本點來更新迴歸係數
    :param data_mat: 輸入數據的數據特徵（除去最後一列）,ndarray
    :param class_labels: 輸入數據的類別標籤（最後一列數據）
    :return: 獲得的最佳迴歸係數
    """
    m, n = np.shape(data_mat)
    alpha = 0.01
    weights = np.ones(n)
    for i in range(m):
        # sum(data_mat[i]*weights)爲了求 f(x)的值， f(x)=a1*x1+b2*x2+..+nn*xn,
        # 此處求出的 h 是一個具體的數值，而不是一個矩陣
        h = sigmoid(sum(data_mat[i] * weights))
        error = class_labels[i] - h
        # 仍是和上面同樣，這個先去看推導，再寫程序
        weights = weights + alpha * error * data_mat[i]
    return weights

In [78]:

def stoc_grad_ascent1(data_mat, class_labels, num_iter=150):
    """
    改進版的隨機梯度上升，使用隨機的一個樣原本更新迴歸係數
    :param data_mat: 輸入數據的數據特徵（除去最後一列）,ndarray
    :param class_labels: 輸入數據的類別標籤（最後一列數據
    :param num_iter: 迭代次數
    :return: 獲得的最佳迴歸係數
    """
    m, n = np.shape(data_mat)
    weights = np.ones(n)
    for j in range(num_iter):
        # 這裏必需要用list，否則後面的del無法使用
        data_index = list(range(m))
        for i in range(m):
            # i和j的不斷增大，致使alpha的值不斷減小，可是不爲0
            alpha = 4 / (1.0 + j + i) + 0.01
            # 隨機產生一個 0～len()之間的一個值
            # random.uniform(x, y) 方法將隨機生成下一個實數，它在[x,y]範圍內,x是這個範圍內的最小值，y是這個範圍內的最大值。
            rand_index = int(np.random.uniform(0, len(data_index)))
            h = sigmoid(np.sum(data_mat[data_index[rand_index]] * weights))
            error = class_labels[data_index[rand_index]] - h
            weights = weights + alpha * error * data_mat[data_index[rand_index]]
            del(data_index[rand_index])
    return weights

In [81]:

weights1 = stoc_grad_ascent1(np.array(data_arr), np.array(label_arr))
weights1

Out[81]:

array([13.97204276,  1.24615972, -1.9102347 ])

In [82]:

plot_best_fit(data_arr, label_arr, weights1)