Python機器學習之Logistic回歸

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

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

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

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

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

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

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

[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行數據（包括特征和標簽），輸出如下。前面我們提到數據中共包含三種鳶尾花（分別標記為setosa、versicolor和virginica），所以這里的標簽 使用的是0，1和2三個數字來分別表示對應的鳶尾花品種，顯然前面5行都屬于標簽為0的鳶尾花。而且一共有150個樣本數據。

[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  

現在我們利用train_test_split函數來對原始數據集進行分類采樣，取其中20%作為測試數據集，取其中80%作為訓練數據集。
[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函數來訓練一個分類器
[python] view plain copy 在CODE上查看代碼片派生到我的代碼片
logreg = linear_model.LogisticRegression(C=1e5, , solver='lbfgs', multi_class='multinomial')  
logreg.fit(X_train, y_train)  
請留意Scikit-Learn文檔中，對于參數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.
參數multi_class的可選值有2個：‘ovr’, ‘multinomial’。多分類問題既可以使用‘ovr’，也可以使用 ‘multinomial’。但是如果你的選項是 ‘ovr’，那么相當于對每個標簽都執行一個二分類處理。Else the loss minimised is the multinomial loss fit acrossthe entire probability distribution. 選項 ‘multinomial’ 則僅適用于將參數solver置為‘lbfgs’時的情況。

然后再利用已經得到的分類器來對測試數據集進行預測
[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))  

預測結果如下，可見總體準確率都在90%以上，分類器執行的還是相當不錯的！
[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()函數和predict()函數來逐條檢視一下Logistic Regression的分類判別結果，請看下面的示例代碼：
[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])  


我們僅僅輸出了前五個測試用例的分類結果，可見這五個樣本的預測結果中前四個都是正確的。
[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  


當然，Logistic Regression的原理網上已有太多資料進行解釋，因此本文的重點顯然并不在于此。但是如果你對該算法的原理比較熟悉，自己實現其中的某些函數也是完全可以的。下面的代碼就演示了筆者自行實現的predict_proba()函數和predict()函數，如果你對此感興趣也不妨試試看。
[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)  

與之前的執行類似，但是這次換成我們自己編寫的函數
[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])  

最后讓我們來對比一下執行結果，可見我們自己實現的函數與直接調用Scikit-Learn中函數所得之結果是完全相同的。
[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  

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

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