Data science & scikit-learn

参考:Documentation of scikit-learn, Machine Learning Map, sklearn user guide



数据导入

数据导入是数据处理、机器学习的第一步。首先对数据进行清洗并存储为csv格式,之后导入数据集为嵌套列表:

import numpy as np
dataset = np.loadtxt('x.csv', delimiter=",")

numpy的更多用法,请查看这里

[TOP]


k近邻

距离表征

Minkowski距离

Minkowski距离也被称为Lp距离,其计算公式为:

\[ d ( {\bf x}, {\bf y} ) = ( \sum_{i=1}^{n} \vert x_i - y_i \vert ^p)^\frac{1}{p} \]

def minkowski_distance(x,y,power):
    if len(x) == len(y):
        return np.power(np.sum(np.power((x-y),power)),(1/(1.0*power)))
    else:
        print "Input should be of equal length"
    return None

这里\( p \geqslant 1 \),当\( p \to \infty \)时,得到Chebyshev距离:

\[ d ( {\bf x}, {\bf y} ) = ( max \vert x_i - y_i \vert ) \]

可以将\( p = 0 \)时的距离看作Hamming距离,Hamming距离最早用于数据传输差错控制编码,它对比两个(相同长度)字符串对应位的取值,将具有不同取值的位的数量定义为Hamming距离,如:

1011101 与 1001001 之间的Hamming距离是 2。

2143896 与 2233796 之间的Hamming距离是 3。

“toned” 与 “roses” 之间的Hamming距离是 3。

对两个字符串进行异或运算,并统计结果为1的个数,这个数就是Hamming距离。

Euclidean距离

令Minkowski距离公式中的p=2,得到Euclidean距离公式:

\[ d ( {\bf x}, {\bf y} ) = ( \sum_{i=1}^{n} \vert x_i - y_i \vert ^2)^\frac{1}{2} \]

def euclidean_distance(x,y):
    return minkowski_distance(x,y,2)

Manhattan距离

令Minkowski距离公式中的p=1,得到Manhattan距离公式:

\[ d ( {\bf x}, {\bf y} ) = \sum_{i=1}^{n} \vert x_i - y_i \vert \]

def euclidean_distance(x,y):
    return minkowski_distance(x,y,1)

cosine距离

将Euclidean距离加1,然后取倒数,就可以得到一个介于0到1之间的值,可以用来表示两个向量相似的程度,1表示完全相同。

与此类似,两个向量夹角的余弦值可以更准确地判断这两个向量的相似程度,它反映的是这两个向量与某一条直线的拟合程度,因此也被称为Pearson相关系数,该系数介于-1到1之间,1表示完全正相关,0表示无关,-1表示完全负相关。Pearson相关系数常常用于文本挖掘,词是轴,文档为向量,Pearson系数就代表了这两个文档之间的相似度。

\[ r ( {\bf x}, {\bf y} ) = \frac{\sum xy - \frac{\sum x \sum y}{n}}{\sqrt{(\sum x^2 - \frac{(\sum x)^2}{n})(\sum y^2 - \frac{(\sum y)^2}{n})}}\]

def cosine_distance(x,y):
    if len(x) == len(y):
        return np.dot(x,y)/np.sqrt(np.dot(x,x)*np.dot(y,y))
    else:
        print "Input should be of equal length"
    return None

对于矩阵a,np.corrcoef(a)np.corrcoef(a,rowvar=0)可以计算各行和各列之间的相关系数,如:

>>> a = [[1, 2, 3, 4, 5],
        [8, 3, 4, 2, 6],
        [9, 5, 7, 2, 5]]
>>> np.corrcoef(a)
array([[ 1.        , -0.32826608, -0.66697297],
      [-0.32826608,  1.        ,  0.80412382],
      [-0.66697297,  0.80412382,  1.        ]])

Jaccard距离

Jaccard距离属于非欧空间距离。两个集合的交集与并集的比值称为Jaccard系数,Jaccard系数可用于衡量聚类算法的性能。1减去Jaccard系数就是Jaccard距离。

def jaccard_distance(x,y):
    set_x = set(x)
    set_y = set(y)
    return 1-len(set_x.intersection(set_y)) / len(set_x.union(set_y))

数据分割

kNN是一种分类算法,其处理的数据集格式通常是这样的:

  1. 每一行是一个实例(或记录、样本等);

  2. 每一列表示一个特征(或属性等);

  3. 最后一列是标签(或类别、标记等)。

对于有监督的kNN算法,需要将导入的数据合理分割为训练集和测试集,为使训练集和测试集中的数据随机分布,首先将数据打乱:

>>> np.random.shuffle(dataset)
>>> train_size = 0.8
>>> test_size = 1-train_size
>>> from sklearn.cross_validation import train_test_split
>>> train, test = train_test_split(dataset,test_size=test_size)

将训练集分解为特征矩阵和标签向量:

for indx in range(len(train[0])):
    X = train[:,:-1]
    y = train[:,-1]

数据的标准化

对于x中不同的特征、属性,取值可能相差很大,如年龄的取值一般为0到100,而行驶里程有可能在0到几万甚至几十万之间取值,这会影响距离计算结果的合理性。因此,有时候,需要对数据集进行归一化处理,也就是进行区间缩放,将所有特征、属性的取值范围映射在0到1之间((X-min)/(max-min)):

from sklearn import preprocessing
X = preprocessing.minmax_scale(X, feature_range=(0, 1))

或者将其取值范围映射到-1到1之间(X/max):

X = preprocessing.maxabs_scale(X)

也可以均值=0,方差=1对数据进行标准化((X-均值)/标准差):

X = preprocessing.scale(X,with_mean=True,with_std=False)

如果设置参数with_std=False,将只使用均值=0进行转换。

注意:上述归一化和标准化都是对特征矩阵中所有列数据分别操作。

还可以对数据进行正则化:

X = preprocessing.normalize_scale(X, norm='l2', return_norm=False)

正则化是将实例缩放到单位范数,也就是使每个实例向量的范数为1,默认为二阶范数l2,也可以选l1max

查找k近邻

指定k值创建kNN分类器,并开始训练:

from sklearn.neighbors import KNeighborsClassifier
knn = KNeighborsClassifier(n_neighbors=10)
knn.fit(x,y)

如果数据量非常大,可以使用partial_fit将训练集分成多次载入内存。

训练完成后将返回knn对象信息:

KNeighborsClassifier(algorithm='auto', leaf_size=30, metric='minkowski',
   metric_params=None, n_jobs=1, n_neighbors=10, p=2,
   weights='uniform')

algorithm, 包含ball_treekd_treebrute以及auto

leaf_size, 整数,用于指定ball_treekd_tree叶节点的数。

metric, 指定距离度量方式。

n_jobs, 指定并行计算线程数量,-1表示CPU的内核数。

n_neighbors指定k值。

p, 指定Minkowski距离中的p值。

weights, 权重,uniform表示等比重,distance表示按距离反比等等。

模型创建完成后就可以用于预测X_test的类别l :

y_predicted = knn.predict(X_test)

计算预测结果的概率和概率的对数:

knn.predict_proba(X_test)
knn.predict_log_proba(X_test)

计算距离测试样本最近的10个点,返回的是这些实例所在的行序号:

knn.kneighbors(X_test,10,False)

第三个参数是return_distance,指定是否返回具体的距离值。

模型性能表征

对照模型预测分类和真实分类,就可以评价该模型的性能:

from sklearn import metrics
print metrics.classification_report(y,y_predicted)

结果是类似这样的报告形式:

        precision    recall  f1-score   support

class 0       0.50      1.00      0.67         1
class 1       0.00      0.00      0.00         1
class 2       1.00      0.67      0.80         3

avg / total       0.70      0.60      0.61         5

查看模型的精度(Accuracy):

knn.score(x,y,sample_weight=None)

查看模型的混淆矩阵(Confusion Matrix):

print metrics.confusion_matrix(y,y_predicted)

混淆矩阵含有丰富的信息:

  Positive Negative  
Positive TP FP P’
Negative FN TN N’
  P N  

可以计算准确率(查准率,Precision):

\[ P = \frac{TP}{TP+FP}=\frac{TP}{P’} \]

print metrics.precision_score(y,y_predicted)

召回率(查全率,Recall):

\[ R = \frac{TP}{TP+FN}=\frac{TP}{P} \]

print metrics.recall_score(y,y_predicted)

精度(准确度,Accuracy):

\[ A = \frac{TP+TN}{P+N} \]

print metrics.accuracy_score(y,y_predicted)

错误率(Error Rate):

\[ E = \frac{FP+FN}{P+N} \]

F1值:

\[ F1 Score = \frac{2PR}{P+R} \]

print metrics.f1_score(y,y_predicted)

此外,metrics模块还有roc_curveauc等方法。

[TOP]


贝叶斯分类器

贝叶斯公式

贝叶斯分类器的理论基础是贝叶斯公式:

\[ P ( c \vert {\bf x } ) = \frac{ P ( {\bf x } \vert c ) P ( c ) }{ P ( {\bf x } ) }\]

这里的P(c)为先验概率,P(x|c)是样本x相对于类标记c的类条件概率,P(x)是用于归一化的因子,对于给定的样本x,P(x)与类标记无关,因此,估计P(c|x)的问题就转化为计算P(x|c)P(c)。贝叶斯分类器与决策树、神经网络、支持向量机等“判别式模型”的不同之处就在于它是一种间接判断模式,被称为“生成式模型”。

如果x中各实例的特征、属性之间相互独立,即每个特征、属性独立影响分类结果,则有

\[ P ( {\bf x } \vert c ) = \prod ^{t}_{i=1} P ( x_i \vert c ) \]

此时的分类器称为朴素贝叶斯分类器。

对\(P( x_i \vert c)\)的不同假设,朴素贝叶斯分类器可分为高斯朴素贝叶斯(Gaussian Naive Bayes),多项式朴素贝叶斯(Multinomial Naive Bayes)和伯努利朴素贝叶斯(Bernoulli Naive Bayes)。

下面尝试使用贝叶斯分类器对文本进行分类过滤。

文本向量化

对文本进行分词(这里以英文为例,中文可使用结巴分词),每个词作为一个元素,组成列表。关于文本处理,可查看这里

sentence = 'today is your last day'
s_vec = [x.strip() for x in sentence.split()]
words_li = [['today', 'is', 'your', 'last', 'day'],
           ['Police', 'have', 'been', 'unable', 'to', 'tell', 'the', 'age'],
           ['We', 'will', 'be', 'keeping', 'this', 'part', 'updated'],
           ['Tis', 'the', 'season', 'to', 'buy', 'votives'],
           ['the', 'lender', 'sends', 'that', 'to', 'the', 'federal', 'bank', 'database'],
           ['Shop', 'service', 'pride', 'items', 'for', 'every', 'branch', 'of', 'service']]

作为训练集,对每一句话的单词列表都进行人工标记,以确定其所属类别:

class_vec = [1,0,0,1,0,1]

其中1表示来自于spam,0是正常文本。

生成一个单词的集合,并转换为列表:

words_set = set([])
for doc in words_li:
    words_set = words_set | set(doc)
all_words_li = list(words_set)

将每一个单词列表都转成向量:

words_matrix = []
for i in range(len(words_li)):
    words_vec = [0]*len(all_words_li)
    for word in words_li[i]:
        if word in all_words_li:
            words_vec[all_words_li.index(word)] = 1
    words_matrix.append(words_vec)

最后将以单词向量为元素,整个训练集就转换为一个嵌套列表(矩阵),记为words_matrix。

高斯朴素贝叶斯

高斯朴素贝叶斯假设\(P( x_i \vert c)\)为正态分布:

\[ P ( x_i \vert c ) = \frac{1}{\sqrt{2 \pi \sigma ^2 _c}} \exp(-\frac{(x_i - \mu _c)^2}{2 \sigma ^2 _c})\]

其中\(\mu _c\)和\(\sigma ^2 _c\)需要从训练集中估计。

from sklearn.naive_bayes import GaussianNB
gnb = GaussianNB()
gnb.fit(words_matrix,class_vec)

对新的输入文本进行分类,如:

new_words_li = ['no', 'news', 'is', 'good', 'news']

将其转为向量new_words_vec,再转为矩阵new_words_matrix后:

gnb.predict(new_words_matrix)

多项式朴素贝叶斯

多项式朴素贝叶斯假设\(P( x_i \vert c)\)为多项式分布:

\[ P ( x_i \vert c ) = \frac{ x_i + \lambda }{ m + n \lambda }\]

from sklearn.naive_bayes import MultinomialNB
mnb = MultinomialNB()
mnb.fit(words_matrix,class_vec)

MultinomialNB()有3个参数,alpha对应于公式中的\(\lambda\),取值通常为1,fit_prior是一个布尔值,表示是否考虑先验概率,当取值为True时,可设定class_prior指定先验概率。

伯努利朴素贝叶斯

伯努利朴素贝叶斯假设\(P( x_i \vert c)\)为二元伯努利分布:

\[ P ( x_i \vert c ) = P ( i \vert c ) x_i + ( 1 - P ( i \vert c ))( 1 - x_i )\]

from sklearn.naive_bayes import BernoulliNB
bnb = BernoulliNB()
bnb.fit(words_matrix,class_vec)

[TOP]


决策树

不纯度表征

信息增益

在给定的样本集合中,每个实例可能属于不同的类别,衡量集合纯度最常用的指标就是熵(entropy),其计算公式为:

\[ Ent ( X ) = - \Sigma ^n _{i=1} P ( x _i ) \log_2 P ( x _i ) \]

熵越小,表示集合的纯度越高,如果集合所有元素属于同一类,则:

\[ Ent ( X ) = - \Sigma ^n _{i=1} 1 \log_2 1 = 0 \]

信息增益可用于表征使用特征(属性)a划分集合时引起的纯度提升,计算公式为:

\[ Gain( X, a ) = Ent ( X ) - \Sigma ^v _{ \upsilon = 1 } \frac{ \vert X^{\upsilon} \vert }{ \vert X \vert } Ent( X^{\upsilon} ) \]

v是特征(属性)a所有可能取值的数目,每一个取值产生一个类。

增益率

增益率(gain ratio)是另外一种衡量集合纯度变化的指标,其计算公式为:

\[ Gain\_Ratio ( X, a ) = \frac{ Gain( X, a )}{ IV(a) } \]

其中:

\[ IV(a) = - \Sigma ^v _{ \upsilon = 1 } \frac{ \vert X^{\upsilon} \vert }{ \vert X \vert } \log_2 \frac{ \vert X^{\upsilon} \vert }{ \vert X \vert } \]

是特征(属性)a的固有熵值(intrinsic value),该值越大,表示特征(属性)a的可能取值数目越多。

基尼指数

基尼值也可用于表征数据集的纯度:

\[ Gini(X) = \Sigma ^{\vert y \vert} _{i=1} \Sigma _{j \neq i} P ( x _i ) P ( x _j ) = 1 - \Sigma ^{\vert y \vert} _{i=1} (P ( x _i ))^2 \]

Gini(X)反映了从数据集中随机抽取两个实例,其类别不一致的概率,Gini(X)越小,数据集的纯度越高。

特征(属性)a的基尼指数的计算公式为:

\[ Gini\_Index (X,a) = \Sigma ^v _{ \upsilon =1} \frac{ \vert X^{\upsilon} \vert }{ \vert X \vert } Gini( X^{\upsilon} ) \]

采用信息增益、增益率和基尼指数的决策树分别为ID3、C4.5和CART算法。下面即将讲述的DecisionTreeClassifier使用的是优化的CART,设置criterion参数为entropy(信息增益),可近似为ID3算法,C4.5基于信息增益率,scikit-learn不能直接实现。

注:ID3中的ID为Iterative Dichotomiser,迭代二分器,CART为Classification and Regression Tree,分类和回归树。

有了上述公式,信息增益、增益率和基尼指数的python实现并不难,对于决策树而言,其可视化更值得关注和研究。

分类决策树

分类决策树所采用的数据集格式与kNN相同,训练集中最后一列为标签,分离特征矩阵和标签向量:

for indx in range(len(train[0])):
    X = train[:,:-1]
    y = train[:,-1]

采用默认的CART算法开始训练:

from sklearn import metrics
from sklearn.tree import DecisionTreeClassifier
dtc = DecisionTreeClassifier()
dtc.fit(X,y)

此时将返回dtc的对象信息:

DecisionTreeClassifier(class_weight=None, criterion='gini', max_depth=None,
    max_features=None, max_leaf_nodes=None,
    min_impurity_decrease=0.0, min_impurity_split=None,
    min_samples_leaf=1, min_samples_split=2,
    min_weight_fraction_leaf=0.0, presort=False, random_state=None,
    splitter='best')

可以看到,可指定的参数非常丰富。

class_weight, 指定各类别实例的权重,防止训练集某些类别的实例过多,导致决策树过于偏向这些类别。默认为’None’;可指定为’balanced’,设定算法自动调整权重;也可以分别指定各个实例的权重,格式为{0:0.9,1:0.1}。

criterion, 可设定为ginientropy,分别对应CART和ID3算法,默认gini

max_depth, 决策树的最大深度,默认为’None’,表示不限制深度,如果实例量或特征数目很多,可设定该值,通常在10-100之间。

max_features, 指定考虑的最大特征数。默认为’None’,表示考虑所有特征;如果为’log2’,表示最多考虑\(log_2N\)个特征(N为数据集特征总数目);如果是’sqrt’或’auto’,表示最多考虑\(\sqrt{N}\)个特征;如果是整数,表示考虑的特征数目;如果是浮点数,表示考虑N乘以该浮点数并取整后的特征数目。如果实例特征数目不多,默认即可;如果特征数目非常多,可设定该参数以控制决策树建模时间。

max_leaf_nodes, 默认为’None’,表示不限制最大的叶子节点数。如果特征数目很多,可以通过交叉验证定一个数值,防止过拟合。

min_impurity_decrease, 浮点数,默认为0,当不纯度的降幅小于该值时终止分类。

min_impurity_split, 默认为2,表示节点实例数少于2时,不会继续尝试选择最优特征进行分类,如果实例数非常大,如10万实例,可考虑设定为10。注:v0.21可能会取消该参数,以min_impurity_decrease替代。

min_samples_leaf, 默认为1,表示叶子节点最少的实例数,可设定为其它整数;如果是浮点数,表示最少实例数占实例总数的百分比。实例量非常大时,可设定该值以提高效率。

min_samples_split, 内部节点再划分所需最小实例数,默认为2。

min_weight_fraction_leaf, 表示叶子节点所有实例权重和的最小值,如果小于该值,会和兄弟节点一起被剪枝;默认是0,表示不考虑权重问题。当实例有较多缺失值或实例分布类别偏差很大时,可考虑设定该值。

presort, 布尔值,默认为’False’,表示不对数据排序,实例量大时,设定为’True’会拖慢建模速度。

splitter, 默认为’best’,适合实例量不大的情况,如果实例量很大,可设定为’random’。

更多参数可查看官方手册。决策树对象还有一些属性和方法,可通过help(dtc)查看。

这时候可以查看当给定x矩阵时,模型预测的标签:

y_predicted = dtc.predict(X_test)

决策树的可视化

决策树经过可视化处理后,其结构会更加清晰,便于理解,决策树可视化一般采用graphviz,需要安装graphviz及相关python模块。

$ sudo apt install graphviz
$ sudo apt install graphviz-dev  #这是和python库链接用的
$ pip install pygraphviz
$ pip install pydotplus  #用于处理dot文件

决策树对象创建完成后可导出为dot文件:

from sklearn import tree
with open('out.dot','w') as f:
    f = tree.export_graphviz(dtc,out_file=f)

dot文件是AutoCAD的图纸文件格式,很多绘图软件都支持,常用于共享数据。在shell中可以快速转为jpg, png, pdf等格式:

$ dot -Tpng out.dot -o out.png

也可以使用pydotplus模块进行处理:

import pydotplus
dot_data = tree.export_graphviz(dtc,out_file=None)
graph = pydotplus.graph_from_dot_data(dot_data)
graph.write_pdf("out.pdf")

模型的保存

决策树可以使用sklearnjoblib模块保存到硬盘或从硬盘导入:

from sklearn.externals import joblib
joblib.dump(dtc,'dtc.model')
dtc = joblib.load('dtc.model')

注:joblib也适用于其他模型、对象的保存和导入。

回归决策树

决策树可以用于回归,此时x为自变量(矩阵),y是因变量(向量),当输入x时返回y,回归就是拟合,目的是寻找y随x的变化规律:

from sklearn import metrics
from sklearn.tree import DecisionTreeRegressor
dtr = DecisionTreeRegressor()
dtr.fit(X,y)

其参数与DecisionTreeClassifier类似,其中,criterion可选msemae,前者为均方差,后者是和均值之差的绝对值之和,默认为mse

完成模型构建后,可以根据给定X_test矩阵,预测相应的y向量:

y_predicted = dtr.predict(X_test)

分类处理的数据集中,y只能是离散值、标签或者’True/False’二元值,而回归处理的数据集中y可以是任意实数。

[TOP]


Logistic回归

Sigmoid函数

对数几率回归虽然带有“回归”两个字,但实际上是一种分类算法。

回归就是根据一系列已知的x和y,去推测函数方程的形式,并计算每一项的系数,拟合出来的直线或曲线应该是连续的,而分类显然不是连续的。

这就要借助于Sigmoid函数(S形函数),最常见的是Logistic函数(对数几率函数),Logistic函数是连续函数,但近似于阶跃函数:

\[ \sigma (z) = \frac{1}{1+e^{-z}} \]

该函数值取值在0到1之间,z=0时,\( \sigma = 0.5\),z减小时,\( \sigma\)快速逼近0,z增大时,\(\sigma\)快速逼近1。

\[ \sigma ‘ (z) = \sigma (z) (1-\sigma(z))\]

可以看到,Logistic函数的导数可以用Logistic函数自身来表示,只要计算出Logistic函数的值,就可以快速得到其导数值。

把回归方程作为z带入Logistic函数中,可以得到一个0到1之间的数值,当该值小于0.5时实例被归入0类,否则归入1类,这样就实现了分类的功能。

范数

范数\(L_p\)是具有“长度”概念的函数,和距离表征有关,向量的范数可以看成是向量的长度。

0阶范数\(L_0\)表示向量到原点的Hamming距离,1阶范数\(L_1\)是向量到原点的Manhattan距离,2阶范数\(L_2\)是Euclidean距离,无穷范数\(L_{\infty}\)是Chebyshev距离。

下面我们将用到\(L_1\)和\(L_2\),其计算公式为:

\[ L_1 = \sum_{i=1}^{n} \vert x_i \vert \]

\[ L_2 = ( \sum_{i=1}^{n} \vert x_i \vert ^2)^\frac{1}{2} \]

如果实例的特征数目很多,容易导致分类模型的过拟合,这时候可以为特征设定权重系数,将不必要的特征变量大幅惩罚降低,\(L_2\)就为这些特征变量分配了近似为0的系数;\(L_1\)则将大量不重要特征的权重系数设置为0,直接丢弃了这些特征变量,使矩阵稀疏化,实现了惩罚的目的。

采用\(L_2\)的Logistic回归称为岭回归,采用\(L_1\)时为LASSO(Least absolute shrinkage and selection operator,最小绝对值缩减和选择)回归。

矩阵求导

A=[a,b,c,d], B=[x,y,z],A对B求导,得到一个(4*3)的矩阵:

\[ \frac{dA}{dB} = \begin{pmatrix} \frac{\partial a}{\partial x}& \frac{\partial b}{\partial x}& \frac{\partial c}{\partial x}& \frac{\partial d}{\partial x}\\
\frac{\partial a}{\partial y}& \frac{\partial b}{\partial y}& \frac{\partial c}{\partial y}& \frac{\partial d}{\partial y}\\
\frac{\partial a}{\partial z}& \frac{\partial b}{\partial z}& \frac{\partial c}{\partial z}& \frac{\partial d}{\partial z} \end{pmatrix} \]

A对B的二阶导数是一个(16*9)的矩阵,称为Hessian矩阵。可见,矩阵求导的计算量和占用内存量都非常大。

梯度下降算法

已知实例的特征矩阵x,标签向量y,确定回归系数\(\omega\)。引入误差\(e\)和步长\(\alpha\),最佳的\(\omega\)需要使\(e^Te\)最小化,为推导方便,引入一个\(\frac{1}{2}\)作为系数,也就是使\(\frac{1}{2}e^Te\)最小化,令:

\[ f(\omega) = \frac{1}{2}e^Te = \frac{1}{2}(x \omega -y)^T(x \omega -y) \]

这里\(f(\omega)\)称为损失函数(loss function)、代价函数(cost function)或残差函数(error function),展开后为

\[ f(\omega) = \frac{1}{2}( \omega ^T x^T x \omega - \omega ^T x ^T y - y^T x \omega + y^Ty) \]

对\(\omega\)求导:

\[ \frac{\partial f (\omega)}{\partial \omega} = - x^T + x^T x \omega = - x^T(y - \omega x) = x^T e \]

可以得到梯度下降算法的迭代式:

\[ \omega := \omega + \frac{\partial f (\omega)}{\partial \omega} = \omega - x^T e \]

引入步长\(\alpha\):

\[ \omega := \omega - \alpha x^T e \]

这里以梯度下降算法阐述了回归系数的确定方法,sklearn的LogisticRegression内置了多种上升或下降算法用于确定迭代方向,如liblinear采用坐标下降算法,lbfgsnewton-cg采用拟牛顿法,sag采用随机平均梯度下降算法,其中,liblinear多用于处理较小的数据集,sag多用于海量数据集。

牛顿法是回归算法中最常用的一种算法,但是由于需要计算二阶导数,生成并存储逆Hessian矩阵,其运算量和需要的内存非常大,因此在工程上用的并不多。BFGS(Broyden-Fletcher-Goldfarb-Shanno)是一种拟牛顿法,只需要计算近似逆Hessian矩阵,但所需内存和计算量仍然很大,sklearn中的lbfgs是Limited-memory BFGS算法,newton-cg是牛顿法的共轭梯度(Conjugate Gradient)近似算法,都属于拟牛顿法,但不需要计算二阶导数,其收敛速度比梯度下降算法快,运算速度比牛顿法快,需要的内存空间比牛顿法小,适用于处理实例数量较大的数据集。

多元分类

从Sigmoid函数的特点可以看到,Logistic回归只能处理二元分类,实际上,也可以处理多元分类,主要有两种方法:ovr即’one vs rest’,把一种类别与其它类别分开,再进一步从其它类别中分出一个类别,经过多次循环就达到了多元分类的目的;multinomial即‘mvm,many vs many’,每次取两个类进行回归,如果共有T类,共需要进行T(T-1)/2次分类,而’ovr’只需要T次,’mvm’算法速度慢一些,但精度更高。

Logistic回归

from sklearn.linear_model import LogisticRegression
lr = LogisticRegression()
lr.fit(X,y)

完成建模后将返回模型的详细参数信息:

LogisticRegression(penalty='l2', dual=False, tol=0.0001, C=1.0,
    fit_intercept=True, intercept_scaling=1, class_weight=None,
    random_state=None, solver='liblinear', max_iter=100,
    multi_class='ovr', verbose=0, warm_start=False, n_jobs=1)

tol, 迭代至误差小于该指定值时终止。

max_iter, 最大迭代次数。

C, 正则化参数,是正则化系数\(\lambda\)的倒数,只能是正数,越小正则越强,正则越强分类越精确,缺点就是容易过拟合。sklearn有LogisticRegressionCV,用交叉验证来选择正则化参数C。

solver, 指定下降算法,包括liblinearlbfgsnewton-cgsag

multi_class, 指定多元分类的实现方法,包括ovrmultinomial

penalty, 指定惩罚基准(正则化参数),包括l1l2

Logistic回归的调参比较复杂,首先可以根据数据量的大小选择算法,如数据量小可优先选择liblinear,数据量非常大可选择sag,数据量较大选lbfgsnewton-cgliblinear不支持’mvm’,如果需要精确分类可考虑其它算法;如果l2过拟合,或数据集的特征很多,如某些生物样本,特征数目多达几千,而实例数目可能只有几百,这时候需要选l1,但l1只有liblinear支持。

[TOP]


支持向量机

核函数

支持向量机(SVM, Support vector machines)是一系列监督学习算法,可用于分类、回归和异常值检测。

支持向量机的推导所涉及的数学知识比较多,过程也比较复杂。用一个超平面将训练集的两类分开,距离该超平面最近的点就是支持向量。为了使分类器通用性、准确性更好,就要寻找一个超平面,使支持向量到该超平面的距离最大,这就是支持向量机的原理。当支持向量确定时,分隔超平面也就确定下来了。

类似于平面上的直线\(y=a x + b\),高维中的超平面可写成\(y= \omega ^T x + b\)的形式,对于SVM分类算法(SVC),两侧的支持向量通过类似于Sigmoid的函数映射到-1到1之间,当映射值小于0时,为-1类,否则为1类,可见,SVC支持向量到分隔超平面的距离总是为1。

低维的不可线性分隔的数据集,可通过核函数(kernel function)”还原”到高维空间中,转换为可线性分隔的问题,这种方法称为核方法。核方法使分类更准确,但是会使计算量增大,工程上引入了松弛变量(slack variable),允许一些数据点处于分隔面的错误一侧,实际操作时,松弛变量通过可惩罚系数(C)来调节,C越小则容错越少。

核函数的选择在很大程度上决定了SVM的性能,常用的核函数有线性核(linear)、多项式核(poly)、高斯核(rbf)、Sigmoid核(sigmoid),这些核函数都包含有一些可调参数,一般需要通过交叉验证来确定。在处理文本数据时可选用线性核,情况不明时先尝试高斯核。

详细的支持向量机原理,可查看这里

支持向量机

sklearn中的svm包括SVCNuSVCLinearSVC分类算法和SVRNuSVRLinearSVR回归算法。其中LinearSVCLinearSVR只支持线性核函数,无法从低维向高维映射,对线性不可分的数据不可用。

当数据可线性拟合时,直接采用LinearSVCLinearSVR,不需要指定核函数;如果无法确定数据的分布状态,采用SVCSVR,需要选择核函数并调参;如果需要调节对训练集的容错率(以支持向量的百分比表示)则选择NuSVCNuSVR

from sklearn import svm
sc = svm.SVC()
sc.fit(X,y)

将返回模型参数信息:

SVC(C=1.0, cache_size=200, class_weight=None, coef0=0.0,
    decision_function_shape=None, degree=3, gamma='auto', kernel='rbf',
    max_iter=-1, probability=False, random_state=None, shrinking=True,
    tol=0.001, verbose=False)

C,惩罚系数,默认为1,NuSVC没有此参数,但是有一个类似的参数nu用于指定训练集训练时的错误率上限,取值范围为(0,1],默认为0.5。

decision_function_shape, 可选择ovoovrovo(one vs one)和Logistic回归部分提到的mvm类似。LinearSVC没有此参数,有类似的multi_class,可选择ovrcrammer_singer,默认ovr

kernel, 核函数,可指定为linearpolyrbfsigmoid,默认为rbfLinearSVC不含此参数。

coef0degreegamma对应核函数的参数r、d、\(\gamma\),一般都由交叉验证来确定。poly含r、d、\(\gamma\)三个参数,rbf含\(\gamma\),sigmoid含r和\(\gamma\)。

SVM回归算法(SVR)类似于上述分类算法,但也有一些不同之处。

SVR(C=1.0, cache_size=200, coef0=0.0, degree=3, epsilon=0.1, gamma='auto',
    kernel='rbf', max_iter=-1, shrinking=True, tol=0.001, verbose=False)

epsilon, 距离误差,SVM回归并没有将支持向量映射到-1到1,此处的\(\epsilon\)表示支持向量到分隔超平面的最大距离误差,NuSVR没有此参数,用Cnu调节错误率。

loss, 损失函数度量,默认epsilon_insensitive,考虑松弛系数,如果选择squared_epsilon_insensitive则不考虑松弛系数,只有LinearSVR含此参数。

[TOP]


k-means聚类

聚类是一种无监督的学习算法,处理的数据没有标记类别,算法会自动构建类别体系。

k-means算法是最简单的聚类算法,基于期望最大化(EM, Expectation Maximization)的思路。k-means随机确定k个初始点(质心,centroid),将数据集中的每个点与距离最近的质心归为一簇,得到k个簇,然后更新每簇质点为该簇的平均值,这个过程多次循环,直至质心不再更新。

轮廓系数

轮廓系数(Silhouette Coefficient)用于表征聚类结果的质量,对于数据集的一个点i,记\(a _i\)=i到簇内其它点的距离的平均值,\(b _i\)=min(i到其它某一簇内所有点的距离的平均值),点i的轮廓系数就是:

\[ S_i = \frac{b _i - a _i}{max(a_i, b_i)} \]

所有点的轮廓系数的平均值,就是聚类结果总的轮廓系数S。S的取值范围为[-1,1],越接近于1,说明聚类效果越好,如果为负数,表明某两个簇之间有重叠,聚类效果很差。

误差平方和(SSE, Sum of Squared Error, 误差向量的二阶范数的平方)也可用于表征聚类结果的性能,SSE越小表明聚类效果越好。

k-means聚类

from sklearn.cluster import KMeans
kmc=KMeans(n_clusters=3)
kmc.fit(X)

返回模型的参数信息:

KMeans(algorithm='auto', copy_x=True, init='k-means++', max_iter=300,
    n_clusters=3, n_init=10, n_jobs=1, precompute_distances='auto',
    random_state=None, tol=0.0001, verbose=0)

algorithm, 包括autofullelkan算法,elkan的效率更高,但不支持离散数据,默认为auto

init, 初始化质心的方法,包括k-means++random或一个ndarray,默认为k-means++

查看聚类结果

kmc.labels_

每一个实例所属的簇都会被标记出来:

array([0, 0, 0, 2, 2, 0, 0, 0, 1, 0, 1], dtype=int32)

查看轮廓系数:

from sklearn.metrics import silhouette_score
silhouette_score(X,kmc.labels_)

[TOP]


bagging

bagging(booststraping aggregating, 自举汇聚),是一种组合算法,即集成学习方法(ensemble learning),bagging每次可放回地抽取训练集的一部分构建基分类器,通过对多个基分类器的结果进行投票决定最终结果。这些基分类器可以是同样的算法,也可以不同,但必须是同类的。bagging是多个基分类器的“并联”。

bagging要求基分类器是不稳定的,也就是数据集微小的变动能够引起分类结果的显著变动,如神经网络、决策树等。

随机森林

随机森林(random forest)是基于决策树的bagging算法,随机森林不仅为每一棵决策树随机抽取部分实例,同时也会随机选取部分特征。随机森林内部会对结果进行评估,因此没有必要对它进行交叉验证或者用一个独立的测试集对结果进行评估,正是由于有这样的优点,随机森林算法曾一度频频出现在各种数据科学竞赛中。

from sklearn.ensemble import RandomForestClassifier
clf = RandomForestClassifier()
clf.fit(X, y)

返回模型参数信息:

RandomForestClassifier(n_estimators=10, criterion='gini', max_depth=None,
    min_samples_split=2, min_samples_leaf=1, min_weight_fraction_leaf=0.0,
    max_features='auto', max_leaf_nodes=None, min_impurity_decrease=0.0,
    min_impurity_split=None, bootstrap=True, oob_score=False, n_jobs=1,
    random_state=None, verbose=0, warm_start=False, class_weight=None)

n_estimators, 森林中树的数量,默认为10。

bootstrap, 默认为True,即有放回的取样。

oob_score, 指定是否采用袋外实例或特征来评估模型。有放回取样中有部分没有被采样的数据,称为袋外数据(‘Out of Bag’),这些数据不属于训练集。默认为‘False’。

列出特征的重要程度:

clf.feature_importances_

[TOP]


Boosting

Boosting是另一种集成学习方法,和bagging的“并联”不同,boosting是基学习器的“串联”,后一个学习器会根据前一个学习器的结果对数据集进行调整,算法的准确性会逐步提升。

如Adaboost(adaptive boosting, 自适应增强)就会为前一个基学习器分错的实例分配更高的权重系数,加权后的数据集会被用来训练下一个学习器,如此循环,直至达到某个足够小的错误率或指定的迭代次数。Boosting的基学习器倾向于采用弱学习器,也就是其预测结果近似于随机猜测,如决策桩等。

梯度提升树

GBDT(Gradient Boosting Decision Tree)又叫 MART(Multiple Additive Regression Tree),是一种boost算法(但不属于Adaboost),由多棵弱决策树组成,所有树的结论累加起来做最终答案。它在被提出之初就和SVM一起被认为是泛化能力较强的算法。

from sklearn.ensemble import GradientBoostClassifier
clf = GradientBoostClassifier()
clf.fit(X, y)

返回模型的参数信息:

GradientBoostingClassifier(loss='deviance', learning_rate=0.1,
    n_estimators=100, subsample=1.0, criterion='friedman_mse',
    min_samples_split=2, min_samples_leaf=1, min_weight_fraction_leaf=0.0,
    max_depth=3, min_impurity_decrease=0.0, min_impurity_split=None,
    init=None, random_state=None, max_features=None, verbose=0,
    max_leaf_nodes=None, warm_start=False, presort='auto')

loss, 损失函数,默认为对数似然损失函数deviance,对二元和多元分类优化较好;如果选用指数损失函数exponential,相当于Adaboost算法。

learning_rate, 弱学习器的权重缩减系数,也就是步长,默认为0.1。

n_estimators, 弱学习器的个数,默认100。

subsample, 取值范围(0,1],指定从数据集中抽取的比例,默认为1.0,也就是使用全部数据集。

criterion, 默认为friedman_mse(优化的mse),也可选mse(平均平方误差)或mae(平均绝对误差)。

init, 说明初始化时采用的弱学习器。

[TOP]


回归

上述分类算法均可用于回归,如:

决策树回归:

from sklearn import tree
model_DecisionTreeRegressor = tree.DecisionTreeRegressor()

线性回归:

from sklearn import linear_model
model_LinearRegression = linear_model.LinearRegression()

SVM回归:

from sklearn import svm
model_SVR = svm.SVR()

KNN回归:

from sklearn import neighbors
model_KNeighborsRegressor = neighbors.KNeighborsRegressor()

随机森林回归:

from sklearn import ensemble
model_RandomForestRegressor = ensemble.RandomForestRegressor(n_estimators=20)

Adaboost回归:

from sklearn import ensemble
model_AdaBoostRegressor = ensemble.AdaBoostRegressor(n_estimators=50)

GBRT回归:

from sklearn import ensemble
model_GradientBoostingRegressor = ensemble.GradientBoostingRegressor(n_estimators=100)

Bagging回归:

from sklearn.ensemble import BaggingRegressor
model_BaggingRegressor = BaggingRegressor()

ExtraTree极端随机树回归:

from sklearn.tree import ExtraTreeRegressor
model_ExtraTreeRegressor = ExtraTreeRegressor()

查看模型回归结果:

import numpy as np
import matplotlib.pyplot as plt
model.fit(X_train,y_train)
score = model.score(X_test, y_test)
result = model.predict(X_test)
plt.figure()
plt.plot(np.arange(len(result)), y_test,'go-',label='true value')
plt.plot(np.arange(len(result)),result,'ro-',label='predict value')
plt.title('score: %f'%score)
plt.legend()
plt.show()

[TOP]


特征工程

通过以上算法实践可以发现数据集的格式均由特征矩阵X和标签向量y构成,将源数据转换为[X, y]形式的过程就是特征工程。实际上,为了使机器学习算法更加准确和高效,还需要对数据做进一步的处理。

数据预处理

无量纲化

主要包括标准化和区间缩放。

标准化就是以均值=0,方差=1对X中的列分别进行转换:

from sklearn.preprocessing import StandardScaler
StandardScaler().fit_transform(X)

区间缩放是利用两个最值对X中的列分别进行区间缩放:

from sklearn.preprocessing import MinMaxScaler
MinMaxScaler().fit_transform(X)

正则化是将X中的行数据转换为标准正态分布:

from sklearn.preprocessing import Normalizer
Normalizer().fit_transform(X)

二值化

设定一个阈值,将数据转换为1和0:

from sklearn.preprocessing import Binarizer
Binarizer(threshold=3).fit_transform(X)

哑编码

将定性的分类转换为OneHot编码,如第一类为001,第二类为010,第三类为100等:

from sklearn.preprocessing import OneHotEncoder
OneHotEncoder().fit_transform(X.reshape((-1,1)))

缺失值处理

默认为取其当前列有效数值的平均值:

from sklearn.preprocessing import Imputer
Imputer(missing_values=-1, strategy='mean', axis=0).fit_transform(X)

其中missing_values用于指定缺失值的表示形式,默认为NaN;strategy指定填充方式,默认为mean,可选medianmost_frequentaxis指定按列(0)或行(1)计算。

特征选择

方差选择法

分别计算各特征的方差,根据阈值,选择方差大于阈值的特征:

from sklearn.feature_selection import VarianceThreshold
VarianceThreshold(threshold=3).fit_transform(X)

卡方检验

计算特征值与标签的相关性,选择最好的若干个特征:

from sklearn.feature_selection import SelectKBest
from sklearn.feature_selection import chi2
SelectKBest(chi2, k=2).fit_transform(X, y)

递归特征消除

使用一个基模型进行多轮训练,消除若干权值系数的特征,再进行下一轮训练:

from sklearn.feature_selection import RFE
from sklearn.linear_model import LogisticRegression
RFE(estimator=LogisticRegression(), n_features_to_select=2).fit_transform(X, y)

模型辅助筛选

使用模型筛选特征:

from sklearn.feature_selection import SelectFromModel
from sklearn.ensemble import GradientBoostingClassifier
SelectFromModel(GradientBoostingClassifier()).fit_transform(X, y)

降维

经过特征选择后,如果特征矩阵仍然特别大,可以进行降维,也就是将数据映射到低维空间。

主成分分析

主成分分析(PCA, )目的是使样本具有最大的发散性:

from sklearn.decomposition import PCA
PCA(n_components=2).fit_transform(X)

n_components为主成分数目。fit_transform不需要标签向量y,可见PCA是一种无监督的算法。

线性判别分析

线性判别分析(LDA, )目标是使样本具有最好的分类特性:

from sklearn.lda import LDA
LDA(n_components=2).fit_transform(X, y)

特征工程需要人的参与,更多地依赖个人的经验,能否尽可能减少人的参与和监督,直接将数据交给机器进行处理,然后自动对不同特征赋予不同权重呢?神经网络在一定程度上可以做到这一点。关于神经网络,请查看这里

文中如有未展开说明的模块,可查看sklearn官方文档

[TOP]