原創：logistic regression實戰（一）：SGD Without lasso

　　logistic regression是分類算法中非常重要的算法，也是非常基礎的算法。logistic regression從整體上考慮樣本預測的精度，用判別學習模型的條件似然進行參數估計，假設樣本遵循iid,參數估計時保證每個樣本的預測值接近真實值的概率最大化。這樣的結果，只能是犧牲一部分的精度來換取另一部分的精度。而svm從局部出發，假設有一個分類平面，找出所有距離分類平面的最近的點(support vector,數量很少)，讓這些點到平面的距離最大化，那么這個分類平面就是最佳分類平面。從這個角度來看待兩個算法，可以得出logistic regression的精度肯定要低於后者。今天主要寫logistic regression的Python代碼。logistic regression的推導過程比較簡單：

　　第一個公式是條件似然函數估計，意思是指定未知常量theta（;表示頻率學派），對於每個輸入feature vector x(i),產生y(i)的概率都最大，取對數是為了求導方便。第二個公式是sigmoid函數的導數，在這里推導出具體的導函數(推導過程非常簡單，復合函數求導法則)，第三個公式是求出的梯度(實際為偏導數組成的向量，梯度與方向導數同方向時取得最大值，相反時取得最小值)。公式3的線性意義為：對於容量為m的樣本(矩陣mxn),權重為nx1的列向量，每一行的樣本數據與權重向量相乘求得預測值，m行樣本的預測值組成mx1的列向量(其數值為exp（-inX）)，其實就是mxn的樣本矩陣左乘權重向量nx1，因為的矩陣的本質是線性變換，相當於每一行樣本數據投影到權重向量。得到預測值的列向量后，與真實值向量(需要轉置）做差值得到新的列向量。最后樣本矩陣倒置，然后左乘這個向量得到nx1的權重向量(梯度)。這個過程是計算批量梯度。當似然函數取得最大時，就是損失函數最小，所以二者是相反的關系，最后結果也就是梯度上升法。

　　對於邏輯回歸的損失函數，也可以從另外一個角度來理解。上面的公式我們看到損失函數和最大似然函數是相反的關系，一般情況下，logistic regression的loss function 可以采用交叉熵的形式，然后取mean。

　　SGD(stochastic gradient descent)在計算的過程中，每次迭代不用計算所用樣本，而是隨機選取一個樣本進行梯度上升(下降)，在工程實踐中可以滿足精度要求，而且時間復雜度比批量要低。對於SGD而言，學習速率alpha的選取，隨着循環次數增加，剛開始的時候，梯度下降應該比較快，到后期的時候，可能會出現在某一個值徘徊的情況，而且下降速率會越來越慢。所以選取一個比較合理的study rate,應該是先選取到的樣本study rate相對較高，后面的相對較小，這樣比較合理。

　　在訓練出模型后，預測時應該使用"五折交叉驗證理論"尋找出最優模型出來(調節參數的過程)，這個過程一般用RMSE（均方根誤差）來衡量，並且可以定義一個基准均方根誤差。在用最優模型預測時取得的RMSE與基准RMSE對比，分析數據結果。

　　本文主要探討logistic regression的SGD算法，並且不考慮正則化。事實上，在高維度的情況下，泛化能力或者正則化技術非常關鍵。在90年代有學者提出lasso后，至今有很多方法實踐了lasso,比如LARS，CD……。下一篇博客將探討lasso技術，並且動手實踐CD算法。接下來，上傳最近寫的SGD Python代碼,首先是引入模塊：logisticRegression.py，這里面定義了兩個class:LogisticRegressionWithSGD,LRModel,還有全局函數RMSE，loadDataSet和sigmoid函數。后面是測試代碼，主要是參數調優。

logisticRegression.py：

'''
Created on 2017/03/15
Logistic Regression without lasso
@author: XueQiang Tong
'''
'''
this algorithm include SGD，batch gradient descent except lasso regularization，So the generalization
ability is relatively weak, follow-up and then write a CD algorithm for lasso.
'''
from numpy import *
import matplotlib.pyplot as plt
import numpy as np;

# load dataset into ndarray(numpy)
def loadDataSet(filepath, seperator="\t"):
    with open(filepath) as fr:
        lines = fr.readlines();
        num_samples = len(lines);
        dimension = len(lines[0].strip().split(seperator));
        dataMat = np.zeros((num_samples, dimension));
        labelMat = [];

    index = 0;
    for line in lines:
        sample = line.strip().split();
        feature = list(map(np.float32, sample[:-1]));
        dataMat[index, 1:] = feature;
        dataMat[index, 0] = 1.0;
        labelMat.append(float(sample[-1]));
        index += 1;

    return dataMat, array(labelMat);

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

# compute rmse
def RMSE(true,predict):
    true_predict = zip(true,predict);
    sub = [];
    for r,p in true_predict:
        sub.append(math.pow((r-p),2));
    return math.sqrt(mean(sub));

def maxstep(dataMat):
    return 2 / abs(np.linalg.det(np.dot(dataMat.transpose(),dataMat)));

class LRModel:
    def __init__(self,weights = np.empty(10),alpha = 0.01,iter = 150):
        self.weights = weights;
        self.alpha = alpha;
        self.iter = iter;

    #predict
    def predict(self,inX):
        prob = sigmoid(dot(inX,self.weights))
        if prob > 0.5: return 1.0
        else: return 0.0

    def getWeights(self):
        return self.weights;

    def __getattr__(self, item):
        if item == 'weights':
            return self.weights;

class LogisticRegressionWithSGD:
    # batch gradient descent
    @classmethod
    def batchGradDescent(cls,data,maxCycles = 500,alpha = 0.001):
        dataMatrix = mat(data[0])             # convert to NumPy matrix
        labelMat = mat(data[1]).transpose() # convert to NumPy matrix
        m,n = shape(dataMatrix)
        weights = ones((n,1))
        for k in range(maxCycles):
            h = sigmoid(dataMatrix*weights)    # compute predict_value
            error = (labelMat - h)              # compute deviation
            weights = weights + alpha * dataMatrix.transpose()* error # update weight

        model = LRModel(weights = weights,alpha = alpha,iter = maxCycles);
        return model;

    '''
    # draw samples of the scatter plot to view the distribution of sample points
    @classmethod
    def viewScatterPlot(cls,weights,dataMat,labelMat):
        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',label='1')
        ax.scatter(xcord2, ycord2, s=30, c='green',label='0')
        plt.legend();
        x = arange(-3.0, 3.0, 0.1)
        y = (-weights[0]-weights[1]*x)/weights[2]
        ax.plot(mat(x), mat(y))
        plt.xlabel('X1'); plt.ylabel('X2');
        plt.show()
    '''
    @classmethod
    def train_nonrandom(cls,data,alpha = 0.01):
        dataMatrix = data[0];
        classLabels = data[1];
        m,n = shape(dataMatrix);
        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];

        model = LRModel(weights = weights,alpha = alpha);
        return model;

    @classmethod
    def train_random(cls,data, numIter=150):
        dataMatrix = data[0];
        classLabels = data[1];
        m,n = shape(dataMatrix);
        weights = ones(n);   #initialize to all ones
        indices = np.arange(len(dataMatrix))
        for iter in range(numIter):
            np.random.shuffle(indices)
            for index in indices:
                alpha = 4 / (1.0 + iter + index) + 0.0001    #apha decreases with iteration, does not
                h = sigmoid(sum(dataMatrix[index]*weights));
                error = classLabels[index] - h;
                weights = weights + alpha * error * array(dataMatrix[index]);

        model = LRModel(weights = weights,iter = numIter);
        return model;

    def colicTest(self):
        data = loadDataSet('G:\\testSet.txt');
        trainWeights = self.stoGradDescent_random(data);
        dataMat = data[0];
        labels = data[1];
        errnums = 0;
        for index in range(dataMat.shape[0]):
            preVal = self.predict(dataMat[index,:],trainWeights);
            if(preVal != labels[index]):
                errnums += 1;
        print('error rate:%.2f' % (errnums/dataMat.shape[0]));
        return errnums;

    def multiTest(self):
        numTests = 10; errorSum=0.0
        for k in range(numTests):
            errorSum += self.colicTest()
        print("after %d iterations the average error rate is: %.2f" % (numTests, errorSum/numTests))

 

 測試代碼，把最優模型保存在npy文件里，以后使用的時候，直接取出來，不用再訓練了。

from logisticRegression import *;
from numpy import *;
import math;

def main():
    dataMat, labels = loadDataSet('G:\\testSet.txt');
    num_samples = dataMat.shape[0];

    num_trains = int(num_samples * 0.6);
    num_validations = int(num_samples * 0.2);
    num_tests = int(num_samples * 0.2);

    data_trains = dataMat[:num_trains, :];
    data_validations = dataMat[num_trains:(num_trains + num_validations), :];
    data_tests = dataMat[(num_trains + num_validations):, :];

    label_trains = labels[:num_trains];
    label_validations = labels[num_trains:(num_trains + num_validations)];
    label_tests = labels[(num_trains + num_validations):];
    '''
    minrmse = (1 << 31) - 1;
    bestModel = LRModel();
    iterList = [10, 20, 30, 80];

    for iter in iterList:
        model = LogisticRegressionWithSGD.train_random((data_trains, label_trains), numIter=iter);
        preVals = zeros(num_validations);
        for i in range(num_validations):
            preVals[i] = model.predict(data_validations[i, :]);

        rmse = RMSE(label_validations, preVals);
        if rmse < minrmse:
            minrmse = rmse;
            bestModel = model;

    print(bestModel.iter, bestModel.weights, minrmse);
    save('D:\\Python\\models\\weights.npy',bestModel.weights);'''

    #用最佳模型預測測試集的評分，並計算和實際評分之間的均方根誤差
    weights = load('D:\\Python\\models\\weights.npy');

    LogisticRegressionWithSGD.viewScatterPlot(weights,dataMat,labels);#顯示散點圖
    model = LRModel();
    model.weights = weights;

    preVals = zeros(num_tests);
    for i in range(num_tests):
        preVals[i] = model.predict(data_tests[i,:]);


    testRMSE = RMSE(label_tests,preVals); #預測產生的均方根誤差
    #用基准偏差衡量最佳模型在測試數據上的預測精度
    tavMean = mean(hstack((label_trains,label_validations)));
    baseRMSE = math.sqrt(mean((label_tests - tavMean) ** 2)) #基准均方根誤差

    improvement = abs(testRMSE - baseRMSE) / baseRMSE * 100;

    print("The best model improves the base line by %% %1.2f" % (improvement));

if __name__ == '__main__':
    main();

 

 

運行結果：
樣本點的散點圖：

最佳迭代次數，權重以及RMSE： 10 [ 12.08509707   1.4723024   -1.86595103] 0.0
The best model improves the base line by % 55.07

另外，由於訓練過程中是隨機選取樣本點，所以迭代次數相同的情況下，權重以及RMSE有可能不同，我們要的是RMSE最小的模型！
算法比較簡單，但是用Python的nump庫實施的時候，有很多注意的細節，只有經過自己仔細的理論推導然后再代碼實施后，才能算基本掌握了一個算法。寫代碼及優化的過程是很費時的，后續還要改進算法，深化理論研究，並且堅持理論與編程結合，切不可眼高手低！