logistic回歸

回歸就是對已知公式的未知參數進行估計。比如已知公式是$y = a*x + b$，未知參數是a和b，利用多真實的(x,y)訓練數據對a和b的取值去自動估計。估計的方法是在給定訓練樣本點和已知的公式后，對於一個或多個未知參數，機器會自動枚舉參數的所有可能取值，直到找到那個最符合樣本點分布的參數（或參數組合）。

logistic分布

設X是連續隨機變量，X服從logistic分布是指X具有下列分布函數和密度函數：

$$F(x)=P(x \le x)=\frac 1 {1+e^{{-(x-\mu)/\gamma}}\\ f(x)=F}{'}(x)=\frac {e^{-(x-\mu)/\gamma}} {\gamma(1+e^{{-(x-\mu)/\gamma})}2}$$

其中，$\mu$為位置參數，$\gamma$為形狀參數。

$f(x)$與$F(x)$圖像如下，其中分布函數是以$(\mu, \frac 1 2)$為中心對陣，$\gamma$越小曲線變化越快。

logistic回歸模型

二項logistic回歸模型如下：

$$P(Y=1|x)=\frac {exp(w \cdot x + b)} {1 + exp(w \cdot x + b)} \
P(Y=0|x)=\frac {1} {1 + exp(w \cdot x + b)}$$

其中，$x \in R^n$是輸入，$Y \in {0,1}$是輸出，w稱為權值向量，b稱為偏置，$w \cdot x$為w和x的內積。

參數估計

假設：

$$P(Y=1|x)=\pi (x), \quad P(Y=0|x)=1-\pi (x)$$

則似然函數為:

$$\prod_{i=1}^N [\pi (x_i)]^{y_i} [1 - \pi(x_i)]^{1-y_i}
$$

求對數似然函數：

$$L(w) = \sum_{i=1}^N [y_i \log{\pi(x_i)} + (1-y_i) \log{(1 - \pi(x_i)})]\
=\sum_{i=1}^N [y_i \log{\frac {\pi (x_i)} {1 - \pi(x_i)}} + \log{(1 - \pi(x_i)})]$$

從而對$L(w)$求極大值，得到$w$的估計值。

求極值的方法可以是梯度下降法，梯度上升法等。

梯度上升確定回歸系數

logistic回歸的sigmoid函數：

$$\sigma (z) = \frac 1 {1 + e^{-z}}$$

假設logstic的函數為:

$$y = w_0 + w_1 x_1 + w_2 x_2 + ... + w_n x_n$$

可簡寫為:

$$y = w_0 + w^T x$$

梯度上升算法是按照上升最快的方向不斷移動，每次都增加$\alpha \nabla_w f(w)$，

$$w = w + \alpha \nabla_w f(w) $$

其中，$\nabla_w f(w)$為函數導數，$\alpha$為增長的步長。

訓練算法

本算法的主要思想根據logistic回歸的sigmoid函數來將函數值映射到有限的空間內，sigmoid函數的取值范圍是0~1，從而可以把數據按照0和1分為兩類。在算法中，首先要初始化所有的w權值為1，每次計算一次誤差並根據誤差調整w權值的大小。

• 每個回歸系數都初始化為1
• 重復N次
  • 計算整個數據集合的梯度
  • 使用$\alpha \cdot \nabla f(x)$來更新w向量
  • 返回回歸系數

#!/usr/bin/env python
# encoding:utf-8

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


def sigmoid(x):
    return 1.0 / (1 + numpy.exp(-x))


def loadData():
    dataMat = []
    laberMat = []
    with open("test.txt", 'r') as f:
        for line in f.readlines():
            arry = line.strip().split()
            dataMat.append([1.0, float(arry[0]), float(arry[1])])
            laberMat.append(float(arry[2]))
    return numpy.mat(dataMat), numpy.mat(laberMat).transpose()


def gradAscent(dataMat, laberMat, alpha=0.001, maxCycle=500):
    """general gradscent"""
    start_time = time.time()
    m, n = numpy.shape(dataMat)
    weights = numpy.ones((n, 1)) 
    for i in range(maxCycle):
        h = sigmoid(dataMat * weights)
        error = laberMat - h 
        weights += alpha * dataMat.transpose() * error
    duration = time.time() - start_time
    print "duration of time:", duration
    return weights


def stocGradAscent(dataMat, laberMat, alpha=0.01):
    start_time = time.time()
    m, n = numpy.shape(dataMat)
    weights = numpy.ones((n, 1)) 
    for i in range(m):
        h = sigmoid(dataMat[i] * weights)
        error = laberMat[i] - h 
        weights += alpha * dataMat[i].transpose() * error
    duration = time.time() - start_time
    print "duration of time:", duration
    return weights


def betterStocGradAscent(dataMat, laberMat, alpha=0.01, numIter=150):
    """better one, use a dynamic alpha"""
    start_time = time.time()
    m, n = numpy.shape(dataMat)
    weights = numpy.ones((n, 1)) 
    for j in range(numIter):
        for i in range(m):
            alpha = 4 / (1 + j + i) + 0.01
            h = sigmoid(dataMat[i] * weights)
            error = laberMat[i] - h 
            weights += alpha * dataMat[i].transpose() * error
    duration = time.time() - start_time
    print "duration of time:", duration
    return weights
    start_time = time.time()


def show(dataMat, laberMat, weights):
    m, n = numpy.shape(dataMat) 
    min_x = min(dataMat[:, 1])[0, 0]
    max_x = max(dataMat[:, 1])[0, 0]
    xcoord1 = []; ycoord1 = []
    xcoord2 = []; ycoord2 = []
    for i in range(m):
        if int(laberMat[i, 0]) == 0:
            xcoord1.append(dataMat[i, 1]); ycoord1.append(dataMat[i, 2]) 
        elif int(laberMat[i, 0]) == 1:
            xcoord2.append(dataMat[i, 1]); ycoord2.append(dataMat[i, 2]) 
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(xcoord1, ycoord1, s=30, c="red", marker="s")
    ax.scatter(xcoord2, ycoord2, s=30, c="green")
    x = numpy.arange(min_x, max_x, 0.1)
    y = (-weights[0] - weights[1]*x) / weights[2]
    ax.plot(x, y)
    plt.xlabel("x1"); plt.ylabel("x2")
    plt.show()

    
if __name__ == "__main__":
    dataMat, laberMat = loadData()
    #weights = gradAscent(dataMat, laberMat, maxCycle=500)
    #weights = stocGradAscent(dataMat, laberMat)
    weights = betterStocGradAscent(dataMat, laberMat, numIter=80)
    show(dataMat, laberMat, weights)

未優化的程序結果如下，

隨機梯度上升算法（降低了迭代的次數，算法較快，但結果不夠准確）結果如下，

對$\alpha$進行優化，動態調整步長（同樣降低了迭代次數，但是由於代碼采用動態調整alpha，提高了結果的准確性），結果如下