1. 程式人生 > >Python實現線性迴歸和邏輯迴歸演算法

Python實現線性迴歸和邏輯迴歸演算法

本文使用python實現了線性迴歸和邏輯迴歸演算法,並使用邏輯迴歸在實際的資料集上預測疝氣病症病馬的死亡率(當然這裡我們的線性迴歸和邏輯迴歸實現是原生的演算法,並沒有考慮正則化係數問題,後期會將其補充完整)。

一、線性迴歸

1.模型表示

h_{\theta}(x)=\theta_{0}+\theta_{1}x_{1}+\theta_{2}x_{2}+...+\theta_{n}x_{n}

2.損失函式

J(\theta)=\frac{1}{2m}\sum_{i=1}^{m}(h_{\theta}(x^{(i)}-y^{(i)})^{2}

3.梯度下降法求解

repeat until convergence {

       \theta_{j}:=\theta_{j}-\alpha \frac{\partial }{\partial \theta_{j}}J(\theta) \quad for (j=1,2,...,n)  

}

根據上式求解:

\theta_{0}:=\theta_{0}-\alpha \frac{1}{m}\sum_{i=1}^{m}(h_{\theta}(x^{(i)})-y^{(i)})

\theta_{1}:=\theta_{1}-\alpha \frac{1}{m}\sum_{i=1}^{m}(h_{\theta}(x^{(i)})-y^{(i)})\cdot x^{(i)}_{1}

\theta_{2}:=\theta_{2}-\alpha \frac{1}{m}\sum_{i=1}^{m}(h_{\theta}(x^{(i)})-y^{(i)})\cdot x^{(i)}_{2}

...

\theta_{n}:=\theta_{n}-\alpha \frac{1}{m}\sum_{i=1}^{m}(h_{\theta}(x^{(i)})-y^{(i)})\cdot x^{(i)}_{n}

update \theta_{0},\theta_{1}.,...,\theta_{n} simultaneously

In other words:

repeat until converagence:{

        \theta_{j}:=\theta_{j}-\alpha \frac{1}{m}\sum_{i=1}^{m}(h_{\theta}(x^{(i)})-y^{(i)})\cdot x^{(i)}_{j} \quad for(j=0,1,2,...,n)

}

二、線性迴歸python程式碼實現

1.線性迴歸的一般過程

(1)收集資料:採用任意的方式採集獲取到資料

(2)準備資料:比如資料的型別是否為數值型。一般結構化的資料最佳

(3)分析資料:對資料進行分析,特徵數量、資料集大小、特徵的選取等

(4)訓練演算法:通過演算法找到最佳的特徵引數

(5)測試演算法:通過準備好的測試資料集來驗證模型的準確性

2.python程式碼

在程式碼的實現過程中,我們只是簡單的使用了兩個特徵引數(x1 和 x2)來寫的,如果資料集有多個引數,需要修改,寫個迴圈即可。(在以後具體的例項中,會將這一塊補充完整。請稍安勿躁。重點是如何使用梯度下降演算法求解的過程見gradAscent方法,並且還額外實現了改進的梯度下降演算法。另外,這裡面不再用具體的數值去測試演算法,寫測試資料太麻煩了。在邏輯迴歸的python實現中,我們用已有的測試資料來檢查編寫的程式是否有問題。測試用的資料到時候也會貼出來,方便大家使用。)

"""
@Time    : 2018/8/24 20:07
@Author  :  xudong
@File    :  
@Contact :  [email protected]
"""
from numpy import mat, shape, ones, random

#載入資料集
def loadDataSet(filePath):
    dataMat = []; labelMat = [] #分別定義一個數據集list和標籤list
    fr = open(filePath) #開啟檔案輸入流
    for line in fr.readlines(): #處理輸入資料
        lineArr = line.strip().split() #按照空格分割每一行資料集 (x1 x2 x3 ...xn y)
        # 將 1.0(x0) x1 x2...值append到定義好的資料集集list中 等同部落格中的 資料集矩陣X
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
        labelMat.append(int(lineArr[2])) #資料集中的真實值 y 這裡需要根據自己的資料集進行修改索引位置
    return dataMat, labelMat


#梯度下降法求解
def gradAscent(dataMatIn, classLabels):
    '''
    :param dataMatIn: 
    :param classLabels: 
    :return: weights 特徵引數(權重係數) theta的值
    '''
    dataMatrix = mat(dataMatIn); # 將讀取完的資料集X轉化成 矩陣X
    labelMat = mat(classLabels).transpose() #轉置 y
    m,n  = shape(dataMatrix) #m是資料集行數,n是特徵的數量
    alpha = 0.001 #學習速率
    maxCycles = 500 #最大迭代次數
    weights = ones((n, 1)) #初始化特徵引數 theta的值
    for k in range(maxCycles):
        error = dataMatrix * weights - labelMat # h(x) - y的值
        weights = weights - alpha * dataMatrix.transpose() * error # w =: w - a(h - y)x
    return weights

#隨機梯度下降演算法
#一次只使用一個樣本點來更新迴歸係數(特徵引數)
def stocGradAscent0(dataMatrix, classLabels):
    m, n = shape(dataMatrix)
    alpha = 0.001
    weights = ones(n)
    for i in range(m):
        error = dataMatrix[i] * weights - classLabels[i]
        weights = weights - alpha * error * dataMatrix[i]
    return weights

#改進的隨機梯度下降演算法
#隨機選取樣本點和調整學習速率來更新迴歸係數
#這個地方隨機選取樣本點的函式式有點問題,請會的幫忙改成一下
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)))
            error = sum(dataMatrix[randIndex] * weights) - classLabels[randIndex]
            weights = weights - alpha * error * dataMatrix[randIndex]
            del(dataMatrix[randIndex])
    return weights

#main 函式
if __name__ == '__main__':
    trainDataPath = "e:/trainSet.txt";
    dataMat, labelMat = loadDataSet(trainDataPath)
    weights = gradAscent(dataMat, labelMat)

(tips:在上述程式碼中,我們使用的迭代次數來代替收斂準則,你也可以修改收斂準則,當損失函式下降的值小於0.001或者更小即為收斂。)

三、邏輯迴歸

(等寫完NG機器學習總結-(四)邏輯迴歸後再來補充 這裡面的知識,下面只給出邏輯迴歸的程式碼實現,並在後續的更新中給出例項,另外邏輯迴歸和線性迴歸的實現程式碼其實差別不大,多的只是用一層sigmoid函式包裝了線性迴歸的模型函式 h)

四、邏輯迴歸python程式碼實現

1.邏輯迴歸一般過程

(1)收集資料:採用任意的方式採集獲取到資料

(2)準備資料:比如資料的型別是否為數值型。一般結構化的資料最佳

(3)分析資料:對資料進行分析,特徵數量、資料集大小、特徵的選取等

(4)訓練演算法:通過演算法找到最佳的特徵引數

(5)測試演算法:通過準備好的測試資料集來驗證模型的準確性

2.python程式碼實現

"""
@Time    : 2018/8/24 20:07
@Author  :  xudong
@File    :  logistic.py
@Contact :  [email protected]
"""
from numpy.ma import exp
from numpy import mat, shape, ones, random

#載入資料集
def loadDataSet(filePath):
    dataMat = []; labelMat = [] #分別定義一個數據集list和標籤list
    fr = open(filePath) #開啟檔案輸入流
    for line in fr.readlines(): #處理輸入資料
        lineArr = line.strip().split() #按照空格分割每一行資料集 (x1 x2 x3 ...xn y)
        # 將 1.0(x0) x1 x2...值append到定義好的資料集集list中 等同部落格中的 資料集矩陣X
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
        labelMat.append(int(lineArr[2])) #資料集中的真實值 y這裡需要根據自己的資料集修改索引位置
    return dataMat, labelMat

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

#梯度下降法求解
def gradAscent(dataMatIn, classLabels):
    '''
    :param dataMatIn: 
    :param classLabels: 
    :return: weights 特徵引數(權重係數) theta的值
    '''
    dataMatrix = mat(dataMatIn); # 將讀取完的資料集X轉化成 矩陣X
    labelMat = mat(classLabels).transpose() #轉置 y
    m,n  = shape(dataMatrix) #m是資料集行數,n是特徵的數量
    alpha = 0.001 #學習速率
    maxCycles = 500 #最大迭代次數
    weights = ones((n, 1)) #初始化特徵引數 theta的值
    for k in range(maxCycles):
        h = sigmoid(dataMatrix * weights)
        error = (h - labelMat) # h(x) - y的值
        weights = weights - alpha * dataMatrix.transpose() * error # w =: w - a(h - y)x
    return weights

#隨機梯度下降演算法
#一次只使用一個樣本點來更新迴歸係數(特徵引數)
def stocGradAscent0(dataMatrix, classLabels):
    m, n = shape(dataMatrix)
    alpha = 0.001
    weights = ones(n)
    for i in range(m):
        h = sigmoid(dataMatrix[i] * weights)
        error = h - classLabels[i]
        weights = weights - alpha * error * dataMatrix[i]
    return weights

#改進的隨機梯度下降演算法
#隨機選取樣本點和調整學習速率來更新迴歸係數
#這個地方隨機選取樣本點的函式式有點問題,請會的幫忙改成一下
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 = h - classLabels[randIndex]
            weights = weights - alpha * error * dataMatrix[randIndex]
            del(dataMatrix[randIndex])
    return weights


#main 函式
if __name__ == '__main__':
    trainDataPath = "e:/trainSet.txt";
    dataMat, labelMat = loadDataSet(trainDataPath)
    weights = gradAscent(dataMat, labelMat)
    print(weights)

這裡面給出trainSet的測試集,有一百條資料trainSet.txt,然後執行程式碼,打印出weights如下圖。

五、例項從疝氣病症預測病馬的死亡率

本小將使用邏輯迴歸來預測患有疝氣病馬的存活問題(分類問題死亡或存活)。這裡的資料包含了368個樣本和28個特徵。原本的資料來自於(UCI機器學習資料庫)。資料集中有30%的值是缺失的。資料集中有缺失值是個非常棘手的問題,很多研究都致力於解決這個問題。這裡不作過多研究,只給出一點簡單的解決方法:

  • 使用特徵的均值來填補缺失值
  • 使用特殊值來填補缺失值,如-1
  • 忽略有缺失值的樣本資料
  • 使用相似樣本的均值來填補缺失值
  • 使用另外的機器學習演算法預測缺失值

處理完資料,現在資料集中特徵的數量只有20個,原本的資料集中類別標籤有三種,分別是“存活”,“已經死亡”,“已經安樂死”,這裡將兩種死亡合併成一種,另外對於本資料集中所有的缺失值都設定為0,處理完的資料集可以在csdn上下載(其csdn上預設下載要一個分,想改為0的,發現改不了,想要資料集的直接發郵件留言)。本次使用的是隨機梯度下降演算法(程式中的stocGradAscent0方法)求解(另外這裡的隨機梯度可能和我們所說的那個隨機梯度下降演算法不一樣,因為在程式中我們使用的當前樣本的一條記錄去更新迭代迴歸引數,而在程式中stocGradAscent1方法才是我們經常提到的隨機梯度下降演算法,隨機選取樣本來更新迴歸引數直到收斂),因為改進的隨機梯度下降演算法(stocGradAscent1)bug沒調出來,沒跑成功(如果有路過的把這個bug改好不妨留下解決方式,十分感謝,初學python有點不熟練)。其整個程式碼如下:

"""
@Time    : 2018/8/24 20:07
@Author  :  xudong
@File    :  logistic.py
@Contact :  [email protected]
"""

from numpy import *

#載入資料集
def loadDataSet(filePath):
    dataMat = []; labelMat = [] #分別定義一個數據集list和標籤list
    fr = open(filePath) #開啟檔案輸入流
    for line in fr.readlines(): #處理輸入資料
        lineArr = line.strip().split('\t') #按照空格分割每一行資料集 (x1 x2 x3 ...xn y)
        # 將 1.0(x0) x1 x2...值append到定義好的資料集集list中 等同部落格中的 資料集矩陣X
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
        labelMat.append(int(lineArr[2])) #資料集中的真實值 y
    return dataMat, labelMat

def loadDataSet2(filePath):
    dataMat = []; labelMat = []
    fr = open(filePath) #開啟檔案輸入流
    for line in fr.readlines(): #處理輸入資料
        lineArr = line.strip().split() #按照空格分割每一行資料集 (x1 x2 x3 ...xn y)
        # 將 x1 x2...值append到定義好的資料集集list中 等同部落格中的 資料集矩陣X
        n = len(lineArr)
        lineCur = []
        for i in range(n-1):
            lineCur.append(float(lineArr[i]))
        dataMat.append(lineCur)
        labelMat.append(float(lineArr[n-1])) #資料集中的真實值 y
    return dataMat, labelMat

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

#梯度下降法求解
def gradAscent(dataMatIn, classLabels):
    '''
    :param dataMatIn: 
    :param classLabels: 
    :return: weights 特徵引數(權重係數) theta的值
    '''
    dataMatrix = mat(dataMatIn); # 將讀取完的資料集X轉化成 矩陣X
    labelMat = mat(classLabels).transpose() #轉置 y
    m,n  = shape(dataMatrix) #m是資料集行數,n是特徵的數量
    alpha = 0.001 #學習速率
    maxCycles = 500 #最大迭代次數
    weights = ones((n, 1)) #初始化特徵引數 theta的值
    for k in range(maxCycles):
        h = sigmoid(dataMatrix * weights)
        error = (h - labelMat) # h(x) - y的值
        weights = weights - alpha * dataMatrix.transpose() * error
        #w =: w - a(h - y)x
    return weights

#隨機梯度下降演算法
#一次只使用一個樣本點來更新迴歸係數(特徵引數)
def stocGradAscent0(dataMatrix, classLabels):
    m, n = shape(dataMatrix)
    alpha = 0.001
    weights = ones(n)
    for i in range(m):
        h = sigmoid(dataMatrix[i] * weights)
        error = classLabels[i] - h
        weights = weights - alpha * error * dataMatrix[i]
    return weights

#改進的隨機梯度下降演算法
#隨機選取樣本點和調整學習速率來更新迴歸係數
def stocGradAscent1(dataMatrix, classLabels, numiter = 150):
    m, n = shape(dataMatrix)
    weights = ones(n)
    for j in range(numiter):
        dataIndex = range(1, 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

#根據死亡的概率判斷死亡
def classifyVector(intX, weights):
    prob = sigmoid(sum(intX * weights))
    if prob > 0.5 : return 1.0
    else : return 0.0

#讀取病馬資料為矩陣,並利用stocAscent演算法求解
def colicTest(trainFile, testFile):
    frTrain = open(trainFile)
    frTest = open(testFile)
    trainingSet = []; trainingLabel = []
    #讀取訓練集資料並儲存為trainingSet
    for line in frTrain.readlines():
        currLine = line.strip().split('\t')
        lineArr = [] # n = len(currLine)
        for i in range(21): #range(n)
            lineArr.append(float(currLine[i]))
        trainingSet.append(lineArr)
        trainingLabel.append(float(currLine[21])) # currLine(n)
    #開始演算法訓練
    trainWeight = stocGradAscent0(array(trainingSet), trainingLabel)
    errorCount = 0; numTestVec = 0.0
    #讀取測試集並儲存為testSet
    for line in frTest.readlines():
        numTestVec += 1.0; #測試集的數量
        currLine = line.strip().split('\t')
        lineArr = []
        for i in range(21):
            lineArr.append(float(currLine[i]))
        if int(classifyVector(array(lineArr), trainWeight)) != int(currLine[21]):
            errorCount += 1
    #計算分類出錯的錯誤率
    errorRate = (float(errorCount) / numTestVec)
    print("the error rate of this test is : %f" %errorRate)
    return errorRate

#多次試驗測試
def multiTest(trainFile, testFile):
    numTests = 10; errorSum = 0.0
    for k in range(numTests):
        errorSum += colicTest(trainFile, testFile)
    print("after %d iterations the average error rate is: %f" %(numTests, errorSum / float(numTests)))

#main 函式
if __name__ == '__main__':
    trainFile = "f:/horseColicTraining.txt"
    testFile = "f:/horseColicTest.txt"
    multiTest(trainFile, testFile)

執行的結果如下圖:

是的這裡的十次執行結果是一樣的,因為我們用的第二種梯度下降演算法,所以結果是一樣的。如果用改進的隨機梯度下降演算法會不一樣的。