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():