[Python]超簡單版logistic-regression 二元分類器實作及範例

跟logistic奮戰了幾天,終於有點眉目的感覺,趁著腦袋瓜還記著的時候記錄下來
借用以前寫過的PLA簡單實作版來修改
http://terrence.logdown.com/posts/290508-python-simple-perceptron-learning-algorithm-implementations

在PLA分類器實作時候因為是實作簡易版,所以該版本對資料預設是線性可分
而到了logistic regression的,即便資料不是線性可分,仍然可以用一些數學公理找到一個還不錯的解
這邊先簡單對logistic regression的概念跟實作做個筆記
順便測試一下recommonmark的跟其他表示數學的工具
https://github.com/cben/mathdown/wiki/math-in-markdown

在做線性回歸的時候,會使用代表一組資料的模型,並使用了最小平方和找出最能帶表資料群的W
輸出會是一個的加權分數,而有時候分類不一定是那麼死硬
我們會希望他是"Soft"(軟性),於是會希望用一個“機率”模型表達某筆資料是圈或叉的機率有多少
因此我們會使用logistic function去將該加權分數限縮到一個0~1之間的數字代表一個“機率”

在當中一個經典的方式就是使用sigmoid function去當作我們的logistic function

丟進一個參數z,當z是無限大的時候輸出是1,負無限大的則是0,剛好0的時候會輸出0.5,借由這個性質可以用他來模擬機率函數
hypothesis就從變成了
是我們的模型對於x猜中是圈圈的機率則可以是看成猜中是叉叉的機率(基於sigmoid有對稱性)
我們做logistic regression的目標就是讓猜中圈圈跟猜中叉叉的機率“最大”
也就是俗稱的Likelihood,以下借用臺大機器學習基石教材的圖

假設我們有N筆資料,y=1是猜圈圈,y=-1是猜叉叉
我們希望這個值越大越好
當是一堆資料的相乘,這個要做最佳化並不是那麼容易
因此我們會做點手腳,把上面公式取log,根據數學公理我們把一堆東西相乘的式子變成一堆東西相加

而上面式子經過整理之後會變成如下

為了讓式子好看,原本要希望猜對機率最大改成希望猜錯的機率最小
裏面的公式帶入轉換後會變成,,一般稱作Cross Entropy Error
要上面的公式最小化的方法,如果是在線性回歸的時候可以直接去算梯度(gradient)等於0

這邊簡化一下過程,上面的Cross Entropy Error微分後的梯度公式長相如下
http://spark.apache.org/docs/latest/mllib-linear-methods.html

這邊直接借用spark那邊的文件,從上面可以看出,由於Cross Entropy Error的梯度方程式並非一個線性方程式
所以我們不能直接使用梯度(gradient)等於0去求最佳解
這邊得回歸PLA的概念:“修正”
我們會使用梯度下降法(gradient descent)去想辦法最佳化
微積分定理告訴我們,一個方程式在某的點的時候,他往梯度方向會上昇最快,反之往梯度反方向的時候會下降最快
用梯度下降的方式去尋找一個竟可能接近最佳解的
更新的公式如下

是第t輪的
是一個常數,代表要更新多少
是我們的梯度
會需要的理由是因為當更新幅度越大,更新得越快但更容易失準(跑過頭)
越小,雖然更新比較慢,但卻能避免跑過頭
我們用去調整更新的步伐去逐步逼近最低點

基本上只要t設得夠大,幾輪更新之後我們可以相信已經有個還不錯而且接近最佳解的

最後實作如下,用之前PLA的例子,但刻意插入一個無法線性分個的點去做比較
原始碼:https://github.com/del680202/MachineLearning-memo/blob/master/src/lr/logistic-regression.py

2016/11/03

最近重新複習了一遍梯度下降法,才發現一開始的實作有錯誤
梯度回傳的應該是一個向量,我卻回傳了實數
或許整個過程算是誤打誤撞跑出一個看起來像結果的東西
一開始是用-1跟1當作標籤去分類
後來上了史丹佛的教材之後,發現用0,1去分類比較符合logistic的精神
因為sigmoid function輸出是0~1的機率,讓logistic去逼近0或1,這樣比較直覺
更新方是是針對每個維度作偏微分後更新,依照上面的公式
假設要更新第j筆權重,那就計算所有資料對權重j的偏微分之後取平均更新,這用個for迴圈可以很簡單作到

不過其實更好的方式是把for的內容收納到矩陣直接用矩陣運算的方式處理掉,這一版本應該是網路上比較容易找到的版本

同時在做logistic的過程中,我們可以順便記錄每次更新權重之後得到新的誤差
若梯度下降有正常運作的話,每次更新後得到的total cost(誤差總和)應該會越來越少逼近到某個曲線
如下圖

若total cost沒有減少的話,通常可能是實作錯誤或是太大,這時候可以試著調低

#!/usr/bin/env python

# encoding: utf-8


import matplotlib.pyplot as plt
import numpy as np
import math

#用0,1去做分類

dataset = np.array([
((1, -0.4, 0.3), 0),
((1, -0.3, -0.1), 0),
((1, -0.2, 0.4), 0),
((1, -0.1, 0.1), 0),
((1, 0.6, -0.5), 0), #非線性分割點

((1, 0.8, 0.7), 1),
((1, 0.9, -0.5), 1),
((1, 0.7, -0.9), 1),
((1, 0.8, 0.2), 1),
((1, 0.4, -0.6), 1)])

#計算機率函數

def sigmoid(z):
    return 1 / (1 + np.exp(-z))
    
#計算平均梯度

def gradient(dataset, w):
    g = np.zeros(len(w))
    for x,y in dataset:
        x = np.array(x)
        error = sigmoid(w.T.dot(x))
        g += (error - y) * x
    return g / len(dataset)

#計算現在的權重的錯誤有多少

def cost(dataset, w):
    total_cost = 0
    for x,y in dataset:
        x = np.array(x)
        error = sigmoid(w.T.dot(x))
        total_cost += abs(y - error)
    return total_cost

def logistic(dataset): #演算法實作

    w = np.zeros(3) #用0 + 0*x1 + 0*x2當作初始設定 

    limit = 10 #更新十次後停下

    eta = 1 #更新幅度

    costs = [] #紀錄每次更新權重後新的cost是多少

    for i in range(limit):
        current_cost = cost(dataset, w)
        print "current_cost=",current_cost
        costs.append(current_cost)
        w = w - eta * gradient(dataset, w)
        eta *= 0.95 #更新幅度,逐步遞減

    #畫出cost的變化曲線,他應該要是不斷遞減 才是正確

    plt.plot(range(limit), costs)
    plt.show()
    return w
#執行


w = logistic(dataset)
#畫圖


ps = [v[0] for v in dataset]
fig = plt.figure()
ax1 = fig.add_subplot(111)

ax1.scatter([v[1] for v in ps[:5]], [v[2] for v in ps[:5]], s=10, c='b', marker="o", label='O')
ax1.scatter([v[1] for v in ps[5:]], [v[2] for v in ps[5:]], s=10, c='r', marker="x", label='X')
l = np.linspace(-2,2)
a,b = -w[1]/w[2], -w[0]/w[2]
ax1.plot(l, a*l + b, 'b-')
plt.legend(loc='upper left');
plt.show()

其結果如下,儘管不是線性分割的資料,卻也還能有個不錯的結果

Stochastic Gradient Descent (SGD)

談到了梯度下降法,順便試著實作看看Stochastic Gradient Descent (SGD) 隨機梯度下降
參考:https://en.wikipedia.org/wiki/Stochastic_gradient_descent

梯度下降法有個小問題,必需每一輪更新都看過“全部資料”來計算梯度做更新
資料量不多的時候或許還好,但是到了到哪邊都是big data的今天
看過“全部資料”這件事就變得花時間了,因此有了SGD這樣的演算法
每次更新只看隨機“一筆”資料的梯度,只看一筆資料可以讓我們往正確方向接近一點點
相比真正的梯度可能一個是走直線一個則是小碎步歪過去,但是步路走的夠多的話,我們可以選擇“相信”SGD能有個還能接受的解

SGD演算法修正:從更新全部資料的梯度,改成只更新一筆資料的梯度

  • 好處:簡單,運算少,應用於big data or online learning, 不用一大批資料,資料可以一筆一筆來做學習
  • 缺點:會有精度上的缺陷

原始碼:https://github.com/del680202/MachineLearning-memo/blob/master/src/lr/stochastic-logistic-regression.py

def logistic(dataset):
    w = np.zeros(3)
    limit = 200
    eta = 0.1
    costs = []
    for i in range(limit):
        current_cost = cost(dataset, w)
        print "current_cost=",current_cost
        costs.append(current_cost)
        w = w - eta * sgd(dataset, w)
        eta = eta * 0.95
    plt.plot(range(limit), costs)
    plt.show()
    return w

先看原來的logistic實作,這邊把gradient函數改成了sgd更新,並且拉大limit(更新次數)
以求讓sgd的精度可以更接近gradient

完全隨機的挑某筆資料算梯度

def sgd(dataset, w):
    #run sgd randomly

    ind = random.randint(0, len(dataset) - 1)
    x, y = dataset[ind]
    x = np.array(x)
    error = sigmoid(w.T.dot(x))
    g = (error - y) * x
    return g

結果如下

在實作過程中有觀察一個的現象,SGD的每次更新後新的cost並沒有總是下降
但卻也能逼近某個低點,先做個記錄
https://notesonml.wordpress.com/2015/07/06/chapter-15-large-scale-machine-learning/

comments powered by Disqus