Python機器學(xué)習(xí)之Logistic回歸

大數(shù)據(jù)時代，數(shù)據(jù)猶如一座巨大的金礦，等待我們?nèi)グl(fā)掘。而機器學(xué)習(xí)和數(shù)據(jù)挖掘的相關(guān)技術(shù)，無疑就是你挖礦探寶的必備利器！工欲善其事，必先利其器。很多初涉該領(lǐng)域的人，最先困惑的一個問題就是，我該選擇哪種“工具”來進行數(shù)據(jù)挖掘和機器學(xué)習(xí)。我這里的工具主要指的是“語言、系統(tǒng)和平臺”。盡管之于機器學(xué)習(xí)而言，語言和平臺從來都不算是核心問題，但是選擇一個你所熟悉的語言和環(huán)境確實可以令你事半功倍。

現(xiàn)在你的選擇可謂相當(dāng)廣泛，例如Matlab、R和Weka都可以用來進行數(shù)據(jù)挖掘和機器學(xué)習(xí)方面的實踐。其中，Matlab是眾所周知的商業(yè)軟件，而R和Weka都是免費軟件。R是應(yīng)用于統(tǒng)計和數(shù)據(jù)分析的首屈一指的計算機語言和平臺，如果你是擁有數(shù)學(xué)或統(tǒng)計學(xué)相關(guān)專業(yè)背景的人，那么使用R來進行數(shù)據(jù)挖掘就是一個相當(dāng)不錯的選擇。我前面有很多介紹利用R語言進行數(shù)據(jù)挖掘的文章可供參考：
在R中使用支持向量機（SVM）進行數(shù)據(jù)挖掘
機器學(xué)習(xí)中的EM算法詳解及R語言實例
Weka的全名是懷卡托智能分析環(huán)境（Waikato Environment for Knowledge Analysis），是一款免費的，非商業(yè)化的，基于Java環(huán)境下開源的機器學(xué)習(xí)（machine learning）以及數(shù)據(jù)挖掘（data mining）軟件。2005年8月，在第11屆ACM SIGKDD國際會議上，懷卡托大學(xué)的Weka小組榮獲了數(shù)據(jù)挖掘和知識探索領(lǐng)域的最高服務(wù)獎，Weka系統(tǒng)得到了廣泛的認可，被譽為數(shù)據(jù)挖掘和機器學(xué)習(xí) 歷史上的里程碑，是現(xiàn)今最完備的數(shù)據(jù)挖掘工具之一。如果你是一個忠實的Java擁護者，那么使用Weka來進行數(shù)據(jù)挖掘就非常明智。

如果你對R和Weka（或Java）都不是很熟悉，那么我今天將向你推薦和介紹另外一個進行機器學(xué)習(xí)和數(shù)據(jù)挖掘的利器——Python。Python是當(dāng)前非常流行的計算機編程語言，相對C、C++來說，Python的門檻極低，可以很輕松的上手和掌握。More importantly，Python用于為數(shù)眾多，而且相當(dāng)完善的軟件包、工具箱來實現(xiàn)功能上的擴展。這一點與R語言來說非常相似（R的擴展包也可以多到超乎你的想象）。

在Python中進行機器學(xué)習(xí)所需要的軟件包主要是Scikit-Learn。Scikit-learn的基本功能主要被分為六個部分，分類，回歸，聚類，數(shù)據(jù)降維，模型選擇，數(shù)據(jù)預(yù)處理。

作為一個范例，我們今天將演示在Python （版本是3.5.1）中基于Scikit-Learn所提供的函數(shù)來實現(xiàn)Logistic Regression。從名字來看，Logistic 回歸 應(yīng)該屬于一種回歸方法（事實上，它也確實可以被用來進行回歸分析），但實際中，它更多的是被用來作為一種“分類器”（Classifier）。而且，機器學(xué)習(xí)中，分類技術(shù)相比于回歸技術(shù)而言也確實是一個更大陣營。
在下面這個示例中，我們會更多的使用僅屬于Scikit-Learn中的函數(shù)來完成任務(wù)。

下面這個例子中的數(shù)據(jù)源于1936年統(tǒng)計學(xué)領(lǐng)域的一代宗師費希爾發(fā)表的一篇重要論文。彼時他收集了三種鳶尾花（分別標(biāo)記為setosa、versicolor和virginica）的花萼和花瓣數(shù)據(jù)。包括花萼的長度和寬度，以及花瓣的長度和寬度。我們將根據(jù)這四個特征（中的兩個）來建立Logistic Regression模型從而實現(xiàn)對三種鳶尾花的分類判別任務(wù)。

首先我們引入一些必要的頭文件，然后讀入數(shù)據(jù)（注意我們僅僅使用前兩個特征）

[python] view plain copy 在CODE上查看代碼片派生到我的代碼片
import numpy as npy  
from sklearn import linear_model, datasets  
from sklearn.cross_validation import train_test_split  
from sklearn.feature_extraction import DictVectorizer  
from sklearn.metrics import accuracy_score, classification_report  
 
iris = datasets.load_iris()  
X = iris.data[:, :2]  # we only take the first two features.  
Y = iris.target  

作為演示，我們來提取其中的前5行數(shù)據(jù)（包括特征和標(biāo)簽），輸出如下。前面我們提到數(shù)據(jù)中共包含三種鳶尾花（分別標(biāo)記為setosa、versicolor和virginica），所以這里的標(biāo)簽 使用的是0，1和2三個數(shù)字來分別表示對應(yīng)的鳶尾花品種，顯然前面5行都屬于標(biāo)簽為0的鳶尾花。而且一共有150個樣本數(shù)據(jù)。

[python] view plain copy 在CODE上查看代碼片派生到我的代碼片
>>> for n in range(5):  
    print(X[n], Y[5])  
 
      
[ 5.1  3.5] 0  
[ 4.9  3. ] 0  
[ 4.7  3.2] 0  
[ 4.6  3.1] 0  
[ 5.   3.6] 0  
 
>>> len(X)  
150  

現(xiàn)在我們利用train_test_split函數(shù)來對原始數(shù)據(jù)集進行分類采樣，取其中20%作為測試數(shù)據(jù)集，取其中80%作為訓(xùn)練數(shù)據(jù)集。
[python] view plain copy
X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size=0.2, random_state=42)  

然后，我們便可以利用LogisticRegression函數(shù)來訓(xùn)練一個分類器
[python] view plain copy 在CODE上查看代碼片派生到我的代碼片
logreg = linear_model.LogisticRegression(C=1e5, , solver='lbfgs', multi_class='multinomial')  
logreg.fit(X_train, y_train)  
請留意Scikit-Learn文檔中，對于參數(shù)solver和multi_class的說明。其中solver的可選值有4個：‘newton-cg’, ‘lbfgs’, ‘liblinear’, ‘sag’。
For small datasets, ‘liblinear’ is a good choice, whereas ‘sag’ is
faster for large ones.
For multiclass problems, only ‘newton-cg’ and ‘lbfgs’ handle
multinomial loss; ‘sag’ and ‘liblinear’ are limited toone-versus-rest schemes.
參數(shù)multi_class的可選值有2個：‘ovr’, ‘multinomial’。多分類問題既可以使用‘ovr’，也可以使用 ‘multinomial’。但是如果你的選項是 ‘ovr’，那么相當(dāng)于對每個標(biāo)簽都執(zhí)行一個二分類處理。Else the loss minimised is the multinomial loss fit acrossthe entire probability distribution. 選項 ‘multinomial’ 則僅適用于將參數(shù)solver置為‘lbfgs’時的情況。

然后再利用已經(jīng)得到的分類器來對測試數(shù)據(jù)集進行預(yù)測
[python] view plain copy 在CODE上查看代碼片派生到我的代碼片
prediction = logreg.predict(X_test)  
print("accuracy score: ")  
print(accuracy_score(y_test, prediction))  
print(classification_report(y_test, prediction))  

預(yù)測結(jié)果如下，可見總體準(zhǔn)確率都在90%以上，分類器執(zhí)行的還是相當(dāng)不錯的！
[python] view plain copy 在CODE上查看代碼片派生到我的代碼片
accuracy score:   
0.9  
             precision    recall  f1-score   support  
 
          0       1.00      1.00      1.00        10  
          1       0.88      0.78      0.82         9  
          2       0.83      0.91      0.87        11  
 
avg / total       0.90      0.90      0.90        30  


In detail, 我們還可以利用predict_proba()函數(shù)和predict()函數(shù)來逐條檢視一下Logistic Regression的分類判別結(jié)果，請看下面的示例代碼：
[python] view plain copy 在CODE上查看代碼片派生到我的代碼片
logreg_proba = logreg.predict_proba(X_test)  
logreg_pred = logreg.predict(X_test)  
for index in range (5):  
    print(logreg_proba[index])  
    print("Predict label:", logreg_pred[index])  
    print("Correct label:", y_test[index])  


我們僅僅輸出了前五個測試用例的分類結(jié)果，可見這五個樣本的預(yù)測結(jié)果中前四個都是正確的。
[python] view plain copy 在CODE上查看代碼片派生到我的代碼片
[  8.86511110e-26   5.64775369e-01   4.35224631e-01]  
Predict label: 1  
Correct label: 1  
[  9.99999942e-01   3.78533501e-08   2.02808786e-08]  
Predict label: 0  
Correct label: 0  
[  9.92889585e-70   8.98623548e-02   9.10137645e-01]  
Predict label: 2  
Correct label: 2  
[  4.40394856e-21   5.97659713e-01   4.02340287e-01]  
Predict label: 1  
Correct label: 1  
[  5.68223824e-43   2.90652338e-01   7.09347662e-01]  
Predict label: 2  
Correct label: 1  


當(dāng)然，Logistic Regression的原理網(wǎng)上已有太多資料進行解釋，因此本文的重點顯然并不在于此。但是如果你對該算法的原理比較熟悉，自己實現(xiàn)其中的某些函數(shù)也是完全可以的。下面的代碼就演示了筆者自行實現(xiàn)的predict_proba()函數(shù)和predict()函數(shù)，如果你對此感興趣也不妨試試看。
[python] view plain copy 在CODE上查看代碼片派生到我的代碼片
class MyLogisticRegression:  
      
    def __init__(self, weights, constants, labels):  
        self.weights = weights  
        self.constants = constants  
        self.labels = labels  
 
    def predict_proba(self,X):  
        proba_list = []  
        len_label = len(self.labels)  
        for n in X: #.toarray():  
            pb = []  
            count = 0  
            for i in range(len_label):  
                value = npy.exp(npy.dot(n, self.weights[i]) + self.constants[i])  
                count = count + value  
                pb.append(value)  
            proba_list.append([x/count for x in pb])  
        return npy.asarray(proba_list)  
      
    def predict(self,X):  
        proba_list = self.predict_proba(X)  
        predicts = []  
        for n in proba_list.tolist():  
            i = n.index(max(n))  
            predicts.append(self.labels[i])  
        return npy.asarray(predicts)  

與之前的執(zhí)行類似，但是這次換成我們自己編寫的函數(shù)
[python] view plain copy 在CODE上查看代碼片派生到我的代碼片
# Print the result based on my functions   
print('\n')  
my_logreg = MyLogisticRegression(logreg.coef_, logreg.intercept_, logreg.classes_)  
my_logreg_proba = my_logreg.predict_proba(X_test)  
my_logreg_pred = my_logreg.predict(X_test)  
for index in range (5):  
    print(my_logreg_proba[index])  
    print("Predict label:",logreg_pred[index])  
    print("Correct label:", y_test[index])  

最后讓我們來對比一下執(zhí)行結(jié)果，可見我們自己實現(xiàn)的函數(shù)與直接調(diào)用Scikit-Learn中函數(shù)所得之結(jié)果是完全相同的。
[python] view plain copy 在CODE上查看代碼片派生到我的代碼片
[  8.86511110e-26   5.64775369e-01   4.35224631e-01]  
Predict label: 1  
Correct label: 1  
[  9.99999942e-01   3.78533501e-08   2.02808786e-08]  
Predict label: 0  
Correct label: 0  
[  9.92889585e-70   8.98623548e-02   9.10137645e-01]  
Predict label: 2  
Correct label: 2  
[  4.40394856e-21   5.97659713e-01   4.02340287e-01]  
Predict label: 1  
Correct label: 1  
[  5.68223824e-43   2.90652338e-01   7.09347662e-01]  
Predict label: 2  
Correct label: 1  

最后需要補充說明的內(nèi)容是，在我們自己編寫的函數(shù)中存在這一句
[python] view plain copy 在CODE上查看代碼片派生到我的代碼片
for n in X: #.toarray():  

請注意我們注釋掉的內(nèi)容，在本篇文章中，我們所使用的數(shù)據(jù)集屬于是標(biāo)準(zhǔn)數(shù)據(jù)集，并不需要我們做Feature extraction。但是在另外一些時候，例如進行自然語言處理時，我們往往要將特征字典轉(zhuǎn)換成一個大的稀疏矩陣，這時我們再編寫上面的函數(shù)時就要使用下面這句來將稀疏矩陣逐行還原
[python] view plain copy 在CODE上查看代碼片派生到我的代碼片
for n in X.toarray():