月度归档：2018年09月

18.1集成学习概述

集成学习的原理正如盲人摸象这个古代寓言所揭示的道理类似：一群盲人第一次遇到大象，想要通过触觉来了解大象。每个人都摸到大象身体的不同部位。但只摸到不同部分，比如鼻子或一条腿。这些人描述的大象是这样的：“它像一条蛇”，“像一根柱子或一棵树”，等等。这些盲人就好比机器学习模型，每个人都是根据自己的假设，并从自己的角度来理解训练数据的多面性。每个人都得到真相的一部分，但不是全部真相。将他们的观点汇集在一起，你就可以得到对数据更加准确的描述。大象是多个部分的组合，每个盲人说的都不完全准确，但综合起来就成了一个相当准确的观点。
集成学习(ensemble learning)可以说是现在非常火爆的机器学习方法了。目前，集成方法在许多著名的机器学习比赛（如 Netflix、KDD 2009 和 Kaggle 比赛）中能够取得很好的名次。
集成学习本身不是一个单独的机器学习算法，而是通过构建并结合多个机器学习器来完成学习任务。也就是我们常说的“博采众长”。集成学习可以用于分类问题集成，回归问题集成，特征选取集成，异常点检测集成等等，可以说所有的机器学习领域都可以看到集成学习的身影。
集成学习的主要思想：对于一个比较复杂的任务，综合许多人的意见来进行决策往往比一家独大好，正所谓集思广益。其过程如下：

一般来说集成学习可以分为三大类：
①用于减少方差的bagging（方差描述的是预测值作为随机变量的离散程度）
②用于减少偏差的boosting（偏差描述的是预测值和真实值之间的差异，即提高拟合能力）
③用于提升预测效果的stacking

18.1.1 Bagging

Bagging是引导聚合的意思。减少一个估计方差的一种方式就是对多个估计进行平均。
Bagging使用装袋采样来获取数据子集训练基础学习器。通常分类任务使用投票的方式集成，而回归任务通过平均的方式集成。
给定一个大小为n的训练集 D，Bagging算法从中均匀、有放回地选出 m个大小为 n' 的子集Di，作为新的训练集。在这 m个训练集上使用分类、回归等算法，则可得到 m个模型，再通过取平均值、取多数票等方法综合产生预测结果，即可得到Bagging的结果。具体如下图：

对于Bagging需要注意的是，每次训练集可以取全部的特征进行训练，也可以随机选取部分特征训练，例如随机森林就是每次随机选取部分特征。
常用的集成算法模型是随机森林和随机树等。
在随机森林中，每个树模型都是装袋采样训练的。另外，特征也是随机选择的，最后对于训练好的树也是随机选择的。
这种处理的结果是随机森林的偏差增加的很少，而由于弱相关树模型的平均，方差也得以降低，最终得到一个方差小，偏差也小的模型。

18.1.2 boosting

Boosting指的是通过算法集合将弱学习器转换为强学习器。boosting的主要原则是训练一系列的弱学习器，所谓弱学习器是指仅比随机猜测好一点点的模型，例如较小的决策树，训练的方式是利用加权的数据。在训练的早期对于错分数据给予较大的权重。
对于训练好的弱分类器，如果是分类任务按照权重进行投票，而对于回归任务进行加权，然后再进行预测。boosting和bagging的区别在于是对加权后的数据利用弱分类器依次进行训练。
boosting是一族可将弱学习器提升为强学习器的算法，这族算法的工作机制类似：
（1）先从初始训练集训练出一个基学习器；
（2）再根据基学习器的表现对训练样本分布进行调整，使得先前基学习器做错的训练样本在后续受到更多关注；
（3）基于调整后的样本分布来训练下一个基学习器；
（4）重复进行上述步骤，直至基学习器数目达到事先指定的值T，最终将这T个基学习器进行加权结合。具体步骤如下图

如果上面这个图还不太直观，大家可参考以下简单示例：
1、假设我们有如下样本图：

图1
2、第一次分类

图2
第2次分类

图3
在图2中被正确测的点有较小的权重（尺寸较小），而被预测错误的点（+）则有较大的权重（尺寸较大）
第3次分类

图4
在图3中被正确测的点有较小的权重（尺寸较小），而被预测错误的点（-）则有较大的权重（尺寸较大）。

第4次综合以上分类

下面描述的算法是最常用的一种boosting算法，叫做AdaBoost，表示自适应boosting。

AdaBoost算法每一轮都要判断当前基学习器是否满足条件，一旦条件不满足，则当前学习器被抛弃，且学习过程停止。
AdaBoost算法中的个体学习器存在着强依赖关系，应用的是串行生成的序列化方法。每一个基生成器的目标，都是为了最小化损失函数。所以，可以说AdaBoost算法注重减小偏差。
由于属于boosting算法族，采用的是加性模型，对每个基学习器的输出结果加权处理，只会得到一个输出预测结果。所以标准的AdaBoost只适用于二分类任务。基于Boosting思想的除AdaBoost外，还有GBDT、XGBoost等。

18.1.3 Stacking

将训练好的所有基模型对训练基进行预测，第j个基模型对第i个训练样本的预测值（概率值或标签）将作为新的训练集中第i个样本的第j个特征值，最后基于新的训练集进行训练。同理，预测的过程也要先经过所有基模型的预测形成新的测试集，最后再对测试集进行预测。如下图所示。

上图可简化为：

其中Meta-Classifier在实际应用中通常使用单层logistic回归模型。
具体算法为：

18.1.3 .1Stacking中元分类层

为何要Meta-Classifier层？设置该层的目的是啥？其原理是什么？等等或许你还不很清楚，没关系。你看了下面这个说明或许就清楚多了。
让我们假设有三个学生名为LR，SVM，KNN，他们争论一个物理问题，他们对正确的答案可能有不同的看法：

他们认为没有办法相互说服他们的情况，他们通过平均估计他们做民主的事情，这个案例是14.他们使用了最简单的集合形式-AKA模型平均。

他们的老师，DL小姐 - 一位数学老师 - 见证了学生们所拥有的论点并决定提供帮助。她问“问题是什么？”，但是学生们拒绝告诉她（因为他们知道提供所有信息对他们不利，除了他们认为她可能会觉得愚蠢他们在争论这么微不足道的事情）。然而，他们确实告诉她这是一个与物理相关的论点。

在这种情况下，教师无法访问初始数据，因为她不知道问题是什么。然而，她确实非常了解学生 - 他们的优点和缺点，她决定她仍然可以帮助解决这个问题。使用历史信息，了解学生过去的表现，以及她知道SVM喜欢物理并且在这个课程中表现优异的事实（加上她的父亲在青年科学家的物理学院工作），她认为最多适当的答案会更像17。

在这种情况下，教师（DL）是元学习者。她使用其他模型（学生）输出的结果作为输入数据。然后，她将其与历史信息结合起来，了解学生过去的表现，以便更好地估计（并帮助解决冲突）。

然而......物理老师RF先生的意见略有不同。他一直在那里，但他一直等到这一刻才行动！ RF先生最近一直在教授LR私人物理课程，以提高他的成绩（错过DL不知道的事情），他认为LR对最终估计的贡献应该更大。因此他声称正确的答案更像是16！

在这种情况下，RF先生也是一个元学习者，他用不同的逻辑处理历史数据 - 他可以访问比DL小姐更多的来源（或不同的历史信息）。

只有校长GBM做出决定，才能解决此争议！ GBM不知道孩子们说了什么，但他很了解他的老师，他更热衷于信任他的物理老师（RF）。他总结答案更像是16.2。

在这种情况下，校长是元级学习者或元学习者的元学习者，并且通过处理他的老师的历史信息，他仍然可以提供比他们的结果的简单平均值更好的估计。
参考文档：
http://blog.kaggle.com/2017/06/15/stacking-made-easy-an-introduction-to-stacknet-by-competitions-grandmaster-marios-michailidis-kazanova/

18.1.3.2Stacking的几种方法

1) 使用分类器产生的特征输出作为meta-classifier的输入
基本使用方法就是，使用前面分类器产生的特征输出作为最后总的meta-classifier的输入数据，以下为利用stacking的基本使用方法实例。
（1）生成数据

from sklearn import datasets

iris = datasets.load_iris()
X, y = iris.data[:, 1:3], iris.target

from sklearn import datasets

iris = datasets.load_iris()

X, y = iris.data[:, 1:3], iris.target

（2）导入需要的库

from sklearn import model_selection
from sklearn.linear_model import LogisticRegression
from sklearn.neighbors import KNeighborsClassifier
from sklearn.naive_bayes import GaussianNB 
from sklearn.ensemble import RandomForestClassifier
from mlxtend.classifier import StackingClassifier
import numpy as np

from sklearn import model_selection

from sklearn.linear_model import LogisticRegression

from sklearn.neighbors import KNeighborsClassifier

from sklearn.naive_bayes import GaussianNB

from sklearn.ensemble import RandomForestClassifier

from mlxtend.classifier import StackingClassifier

import numpy as np

（3）训练各种基模型

clf1 = KNeighborsClassifier(n_neighbors=1)
clf2 = RandomForestClassifier(random_state=1)
clf3 = GaussianNB()
lr = LogisticRegression()
sclf = StackingClassifier(classifiers=[clf1, clf2, clf3], 
                          meta_classifier=lr)

print('3-fold cross validation:\n')

for clf, label in zip([clf1, clf2, clf3, sclf], 
                      ['KNN', 
                       'Random Forest', 
                       'Naive Bayes',
                       'StackingClassifier']):

    scores = model_selection.cross_val_score(clf, X, y, 
                                              cv=3, scoring='accuracy')
    print("Accuracy: %0.2f (+/- %0.2f) [%s]" 
          % (scores.mean(), scores.std(), label))

clf1 = KNeighborsClassifier(n_neighbors=1)

clf2 = RandomForestClassifier(random_state=1)

clf3 = GaussianNB()

lr = LogisticRegression()

sclf = StackingClassifier(classifiers=[clf1, clf2, clf3],

meta_classifier=lr)

print('3-fold cross validation:\n')

for clf, label in zip([clf1, clf2, clf3, sclf],

['KNN',

'Random Forest',

'Naive Bayes',

'StackingClassifier']):

scores = model_selection.cross_val_score(clf, X, y,

cv=3, scoring='accuracy')

print("Accuracy: %0.2f (+/- %0.2f) [%s]"

% (scores.mean(), scores.std(), label))

运行结果
3-fold cross validation:

Accuracy: 0.91 (+/- 0.01) [KNN]
Accuracy: 0.91 (+/- 0.06) [Random Forest]
Accuracy: 0.92 (+/- 0.03) [Naive Bayes]
Accuracy: 0.95 (+/- 0.03) [StackingClassifier]
（4）可视化结果

%matplotlib inline
import matplotlib.pyplot as plt
from mlxtend.plotting import plot_decision_regions
import matplotlib.gridspec as gridspec
import itertools

gs = gridspec.GridSpec(2, 2)

fig = plt.figure(figsize=(10,8))

for clf, lab, grd in zip([clf1, clf2, clf3, sclf], 
                         ['KNN', 
                          'Random Forest', 
                          'Naive Bayes',
                          'StackingClassifier'],
                          itertools.product([0, 1], repeat=2)):

    clf.fit(X, y)
    ax = plt.subplot(gs[grd[0], grd[1]])
    fig = plot_decision_regions(X=X, y=y, clf=clf)
    plt.title(lab)

%matplotlib inline

import matplotlib.pyplot as plt

from mlxtend.plotting import plot_decision_regions

import matplotlib.gridspec as gridspec

import itertools

gs = gridspec.GridSpec(2, 2)

fig = plt.figure(figsize=(10,8))

for clf, lab, grd in zip([clf1, clf2, clf3, sclf],

['KNN',

'Random Forest',

'Naive Bayes',

'StackingClassifier'],

itertools.product([0, 1], repeat=2)):

clf.fit(X, y)

ax = plt.subplot(gs[grd[0], grd[1]])

fig = plot_decision_regions(X=X, y=y, clf=clf)

plt.title(lab)

运行结果

使用网格方法选择超参数

from sklearn.linear_model import LogisticRegression
from sklearn.neighbors import KNeighborsClassifier
from sklearn.naive_bayes import GaussianNB 
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import GridSearchCV
from mlxtend.classifier import StackingClassifier

# Initializing models

clf1 = KNeighborsClassifier(n_neighbors=1)
clf2 = RandomForestClassifier(random_state=1)
clf3 = GaussianNB()
lr = LogisticRegression()
sclf = StackingClassifier(classifiers=[clf1, clf2, clf3], 
                          meta_classifier=lr)

params = {'kneighborsclassifier__n_neighbors': [1, 5],
          'randomforestclassifier__n_estimators': [10, 50],
          'meta-logisticregression__C': [0.1, 10.0]}

grid = GridSearchCV(estimator=sclf, 
                    param_grid=params, 
                    cv=5,
                    refit=True)
grid.fit(X, y)

cv_keys = ('mean_test_score', 'std_test_score', 'params')

for r, _ in enumerate(grid.cv_results_['mean_test_score']):
    print("%0.3f +/- %0.2f %r"
          % (grid.cv_results_[cv_keys[0]][r],
             grid.cv_results_[cv_keys[1]][r] / 2.0,
             grid.cv_results_[cv_keys[2]][r]))

print('Best parameters: %s' % grid.best_params_)
print('Accuracy: %.2f' % grid.best_score_)

from sklearn.linear_model import LogisticRegression

from sklearn.neighbors import KNeighborsClassifier

from sklearn.naive_bayes import GaussianNB

from sklearn.ensemble import RandomForestClassifier

from sklearn.model_selection import GridSearchCV

from mlxtend.classifier import StackingClassifier

# Initializing models

clf1 = KNeighborsClassifier(n_neighbors=1)

clf2 = RandomForestClassifier(random_state=1)

clf3 = GaussianNB()

lr = LogisticRegression()

sclf = StackingClassifier(classifiers=[clf1, clf2, clf3],

meta_classifier=lr)

params = {'kneighborsclassifier__n_neighbors': [1, 5],

'randomforestclassifier__n_estimators': [10, 50],

'meta-logisticregression__C': [0.1, 10.0]}

grid = GridSearchCV(estimator=sclf,

param_grid=params,

cv=5,

refit=True)

grid.fit(X, y)

cv_keys = ('mean_test_score', 'std_test_score', 'params')

for r, _ in enumerate(grid.cv_results_['mean_test_score']):

print("%0.3f +/- %0.2f %r"

% (grid.cv_results_[cv_keys[0]][r],

grid.cv_results_[cv_keys[1]][r] / 2.0,

grid.cv_results_[cv_keys[2]][r]))

print('Best parameters: %s' % grid.best_params_)

print('Accuracy: %.2f' % grid.best_score_)

运行结果
0.667 +/- 0.00 {'kneighborsclassifier__n_neighbors': 1, 'meta-logisticregression__C': 0.1, 'randomforestclassifier__n_estimators': 10}
0.667 +/- 0.00 {'kneighborsclassifier__n_neighbors': 1, 'meta-logisticregression__C': 0.1, 'randomforestclassifier__n_estimators': 50}
0.927 +/- 0.02 {'kneighborsclassifier__n_neighbors': 1, 'meta-logisticregression__C': 10.0, 'randomforestclassifier__n_estimators': 10}
0.913 +/- 0.03 {'kneighborsclassifier__n_neighbors': 1, 'meta-logisticregression__C': 10.0, 'randomforestclassifier__n_estimators': 50}
0.667 +/- 0.00 {'kneighborsclassifier__n_neighbors': 5, 'meta-logisticregression__C': 0.1, 'randomforestclassifier__n_estimators': 10}
0.667 +/- 0.00 {'kneighborsclassifier__n_neighbors': 5, 'meta-logisticregression__C': 0.1, 'randomforestclassifier__n_estimators': 50}
0.933 +/- 0.02 {'kneighborsclassifier__n_neighbors': 5, 'meta-logisticregression__C': 10.0, 'randomforestclassifier__n_estimators': 10}
0.940 +/- 0.02 {'kneighborsclassifier__n_neighbors': 5, 'meta-logisticregression__C': 10.0, 'randomforestclassifier__n_estimators': 50}
Best parameters: {'kneighborsclassifier__n_neighbors': 5, 'meta-logisticregression__C': 10.0, 'randomforestclassifier__n_estimators': 50}
Accuracy: 0.94

2）使用类别概率值作为meta-classfier的输入
另一种使用第一层基本分类器产生的类别概率值作为meta-classfier的输入，这种情况下需要将StackingClassifier的参数设置为 use_probas=True。如果将参数设置为 average_probas=True，那么这些基分类器对每一个类别产生的概率值会被平均，否则会拼接。
例如有两个基分类器产生的概率输出为：
classifier 1: [0.2, 0.5, 0.3]
classifier 2: [0.3, 0.4, 0.4]
1) average = True :
产生的meta-feature 为：[0.25, 0.45, 0.35]
2) average = False:
产生的meta-feature为：[0.2, 0.5, 0.3, 0.3, 0.4, 0.4]
以下为具体实例

from sklearn import datasets
 
iris = datasets.load_iris()
X, y = iris.data[:, 1:3], iris.target
 
from sklearn import model_selection
from sklearn.linear_model import LogisticRegression
from sklearn.neighbors import KNeighborsClassifier
from sklearn.naive_bayes import GaussianNB 
from sklearn.ensemble import RandomForestClassifier
from mlxtend.classifier import StackingClassifier
import numpy as np
 
clf1 = KNeighborsClassifier(n_neighbors=1)
clf2 = RandomForestClassifier(random_state=1)
clf3 = GaussianNB()
lr = LogisticRegression()
sclf = StackingClassifier(classifiers=[clf1, clf2, clf3],
                          use_probas=True,
                          average_probas=False,
                          meta_classifier=lr)
 
print('3-fold cross validation:\n')
 
for clf, label in zip([clf1, clf2, clf3, sclf], 
                      ['KNN', 
                       'Random Forest', 
                       'Naive Bayes',
                       'StackingClassifier']):
 
    scores = model_selection.cross_val_score(clf, X, y, 
                                              cv=3, scoring='accuracy')
    print("Accuracy: %0.2f (+/- %0.2f) [%s]" 
          % (scores.mean(), scores.std(), label))

from sklearn import datasets

iris = datasets.load_iris()

X, y = iris.data[:, 1:3], iris.target

from sklearn import model_selection

from sklearn.linear_model import LogisticRegression

from sklearn.neighbors import KNeighborsClassifier

from sklearn.naive_bayes import GaussianNB

from sklearn.ensemble import RandomForestClassifier

from mlxtend.classifier import StackingClassifier

import numpy as np

clf1 = KNeighborsClassifier(n_neighbors=1)

clf2 = RandomForestClassifier(random_state=1)

clf3 = GaussianNB()

lr = LogisticRegression()

sclf = StackingClassifier(classifiers=[clf1, clf2, clf3],

use_probas=True,

average_probas=False,

meta_classifier=lr)

print('3-fold cross validation:\n')

for clf, label in zip([clf1, clf2, clf3, sclf],

['KNN',

'Random Forest',

'Naive Bayes',

'StackingClassifier']):

scores = model_selection.cross_val_score(clf, X, y,

cv=3, scoring='accuracy')

print("Accuracy: %0.2f (+/- %0.2f) [%s]"

% (scores.mean(), scores.std(), label))

运行结果
3-fold cross validation:

Accuracy: 0.91 (+/- 0.01) [KNN]
Accuracy: 0.91 (+/- 0.06) [Random Forest]
Accuracy: 0.92 (+/- 0.03) [Naive Bayes]
Accuracy: 0.94 (+/- 0.03) [StackingClassifier]

显然，用stacking方法的精度（Accuracy: 0.94）明显好于单个模型的精度。
3）使用堆叠分类及网格搜索
使用堆叠分类及网格搜索（Stacked Classification and GridSearch）方法，要为scikit-learn网格搜索设置参数网格，我们只需在参数网格中提供估算器的名称 - 在meta-regressor的特殊情况下，我们附加'meta-'前缀即可，以下为代码实例。

from sklearn.linear_model import LogisticRegression
from sklearn.neighbors import KNeighborsClassifier
from sklearn.naive_bayes import GaussianNB 
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import GridSearchCV
from mlxtend.classifier import StackingClassifier

# Initializing models

clf1 = KNeighborsClassifier(n_neighbors=1)
clf2 = RandomForestClassifier(random_state=1)
clf3 = GaussianNB()
lr = LogisticRegression()
sclf = StackingClassifier(classifiers=[clf1, clf2, clf3], 
                          meta_classifier=lr)

params = {'kneighborsclassifier__n_neighbors': [1, 5],
          'randomforestclassifier__n_estimators': [10, 50],
          'meta-logisticregression__C': [0.1, 10.0]}

grid = GridSearchCV(estimator=sclf, 
                    param_grid=params, 
                    cv=5,
                    refit=True)
grid.fit(X, y)

cv_keys = ('mean_test_score', 'std_test_score', 'params')

for r, _ in enumerate(grid.cv_results_['mean_test_score']):
    print("%0.3f +/- %0.2f %r"
          % (grid.cv_results_[cv_keys[0]][r],
             grid.cv_results_[cv_keys[1]][r] / 2.0,
             grid.cv_results_[cv_keys[2]][r]))

print('Best parameters: %s' % grid.best_params_)
print('Accuracy: %.2f' % grid.best_score_)

from sklearn.linear_model import LogisticRegression

from sklearn.neighbors import KNeighborsClassifier

from sklearn.naive_bayes import GaussianNB

from sklearn.ensemble import RandomForestClassifier

from sklearn.model_selection import GridSearchCV

from mlxtend.classifier import StackingClassifier

# Initializing models

clf1 = KNeighborsClassifier(n_neighbors=1)

clf2 = RandomForestClassifier(random_state=1)

clf3 = GaussianNB()

lr = LogisticRegression()

sclf = StackingClassifier(classifiers=[clf1, clf2, clf3],

meta_classifier=lr)

params = {'kneighborsclassifier__n_neighbors': [1, 5],

'randomforestclassifier__n_estimators': [10, 50],

'meta-logisticregression__C': [0.1, 10.0]}

grid = GridSearchCV(estimator=sclf,

param_grid=params,

cv=5,

refit=True)

grid.fit(X, y)

cv_keys = ('mean_test_score', 'std_test_score', 'params')

for r, _ in enumerate(grid.cv_results_['mean_test_score']):

print("%0.3f +/- %0.2f %r"

% (grid.cv_results_[cv_keys[0]][r],

grid.cv_results_[cv_keys[1]][r] / 2.0,

grid.cv_results_[cv_keys[2]][r]))

print('Best parameters: %s' % grid.best_params_)

print('Accuracy: %.2f' % grid.best_score_)

运行结果
0.667 +/- 0.00 {'kneighborsclassifier__n_neighbors': 1, 'meta-logisticregression__C': 0.1, 'randomforestclassifier__n_estimators': 10}
0.667 +/- 0.00 {'kneighborsclassifier__n_neighbors': 1, 'meta-logisticregression__C': 0.1, 'randomforestclassifier__n_estimators': 50}
0.967 +/- 0.01 {'kneighborsclassifier__n_neighbors': 1, 'meta-logisticregression__C': 10.0, 'randomforestclassifier__n_estimators': 10}
0.967 +/- 0.01 {'kneighborsclassifier__n_neighbors': 1, 'meta-logisticregression__C': 10.0, 'randomforestclassifier__n_estimators': 50}
0.667 +/- 0.00 {'kneighborsclassifier__n_neighbors': 5, 'meta-logisticregression__C': 0.1, 'randomforestclassifier__n_estimators': 10}
0.667 +/- 0.00 {'kneighborsclassifier__n_neighbors': 5, 'meta-logisticregression__C': 0.1, 'randomforestclassifier__n_estimators': 50}
0.967 +/- 0.01 {'kneighborsclassifier__n_neighbors': 5, 'meta-logisticregression__C': 10.0, 'randomforestclassifier__n_estimators': 10}
0.967 +/- 0.01 {'kneighborsclassifier__n_neighbors': 5, 'meta-logisticregression__C': 10.0, 'randomforestclassifier__n_estimators': 50}
Best parameters: {'kneighborsclassifier__n_neighbors': 1, 'meta-logisticregression__C': 10.0, 'randomforestclassifier__n_estimators': 10}
Accuracy: 0.97
最一句是各参数最佳匹配模型的结果，显然，这个精度高于其他情况的精度。
如果我们计划多次使用回归算法，我们需要做的就是在参数网格中添加一个额外的数字后缀，如下所示：

from sklearn.model_selection import GridSearchCV

# Initializing models

clf1 = KNeighborsClassifier(n_neighbors=1)
clf2 = RandomForestClassifier(random_state=1)
clf3 = GaussianNB()
lr = LogisticRegression()
sclf = StackingClassifier(classifiers=[clf1, clf1, clf2, clf3], 
                          meta_classifier=lr)

params = {'kneighborsclassifier-1__n_neighbors': [1, 5],
          'kneighborsclassifier-2__n_neighbors': [1, 5],
          'randomforestclassifier__n_estimators': [10, 50],
          'meta-logisticregression__C': [0.1, 10.0]}

grid = GridSearchCV(estimator=sclf, 
                    param_grid=params, 
                    cv=5,
                    refit=True)
grid.fit(X, y)

cv_keys = ('mean_test_score', 'std_test_score', 'params')

for r, _ in enumerate(grid.cv_results_['mean_test_score']):
    print("%0.3f +/- %0.2f %r"
          % (grid.cv_results_[cv_keys[0]][r],
             grid.cv_results_[cv_keys[1]][r] / 2.0,
             grid.cv_results_[cv_keys[2]][r]))

print('Best parameters: %s' % grid.best_params_)
print('Accuracy: %.2f' % grid.best_score_)

from sklearn.model_selection import GridSearchCV

# Initializing models

clf1 = KNeighborsClassifier(n_neighbors=1)

clf2 = RandomForestClassifier(random_state=1)

clf3 = GaussianNB()

lr = LogisticRegression()

sclf = StackingClassifier(classifiers=[clf1, clf1, clf2, clf3],

meta_classifier=lr)

params = {'kneighborsclassifier-1__n_neighbors': [1, 5],

'kneighborsclassifier-2__n_neighbors': [1, 5],

'randomforestclassifier__n_estimators': [10, 50],

'meta-logisticregression__C': [0.1, 10.0]}

grid = GridSearchCV(estimator=sclf,

param_grid=params,

cv=5,

refit=True)

grid.fit(X, y)

cv_keys = ('mean_test_score', 'std_test_score', 'params')

for r, _ in enumerate(grid.cv_results_['mean_test_score']):

print("%0.3f +/- %0.2f %r"

% (grid.cv_results_[cv_keys[0]][r],

grid.cv_results_[cv_keys[1]][r] / 2.0,

grid.cv_results_[cv_keys[2]][r]))

print('Best parameters: %s' % grid.best_params_)

print('Accuracy: %.2f' % grid.best_score_)

运行结果
0.667 +/- 0.00 {'kneighborsclassifier-1__n_neighbors': 1, 'kneighborsclassifier-2__n_neighbors': 1, 'meta-logisticregression__C': 0.1, 'randomforestclassifier__n_estimators': 10}
0.667 +/- 0.00 {'kneighborsclassifier-1__n_neighbors': 1, 'kneighborsclassifier-2__n_neighbors': 1, 'meta-logisticregression__C': 0.1, 'randomforestclassifier__n_estimators': 50}
0.967 +/- 0.01 {'kneighborsclassifier-1__n_neighbors': 1, 'kneighborsclassifier-....
0.667 +/- 0.00 {'kneighborsclassifier-1__n_neighbors': 5, 'kneighborsclassifier-2__n_neighbors': 1, 'meta-logisticregression__C': 0.1, 'randomforestclassifier__n_estimators': 50}
0.967 +/- 0.01 {'kneighborsclassifier-1__n_neighbors': 5, 'kneighborsclassifier-2__n_neighbors': 1, 'meta-logisticregression__C': 10.0, 'randomforestclassifier__n_estimators': 10}
.......................................
0.967 +/- 0.01 {'kneighborsclassifier-1__n_neighbors': 5, 'kneighborsclassifier-2__n_neighbors': 5, 'meta-logisticregression__C': 10.0, 'randomforestclassifier__n_estimators': 50}
Best parameters: {'kneighborsclassifier-1__n_neighbors': 1, 'kneighborsclassifier-2__n_neighbors': 1, 'meta-logisticregression__C': 10.0, 'randomforestclassifier__n_estimators': 10}
Accuracy: 0.97

StackingClassifier还可以对分类器参数进行网格搜索。但是，由于目前scikit-learn中GridSearchCV的实现，不可能同时搜索不同分类器和分类器参数。例如，虽然以下参数字典有效。
4）给不同及分类器不同特征
给不同及分类器不同特征是对训练基中的特征维度进行操作的，这次不是给每一个基分类器全部的特征，而是给不同的基分类器分不同的特征，即比如基分类器1训练前半部分特征，基分类器2训练后半部分特征（可以通过sklearn 的pipelines 实现）。最终通过StackingClassifier组合起来。以下为代码实例。

from sklearn.datasets import load_iris
from mlxtend.classifier import StackingClassifier
from mlxtend.feature_selection import ColumnSelector
from sklearn.pipeline import make_pipeline
from sklearn.linear_model import LogisticRegression

iris = load_iris()
X = iris.data
y = iris.target

pipe1 = make_pipeline(ColumnSelector(cols=(0, 2)),
                      LogisticRegression())
pipe2 = make_pipeline(ColumnSelector(cols=(1, 2, 3)),
                      LogisticRegression())

sclf = StackingClassifier(classifiers=[pipe1, pipe2], 
                          meta_classifier=LogisticRegression())

sclf.fit(X, y)

from sklearn.datasets import load_iris

from mlxtend.classifier import StackingClassifier

from mlxtend.feature_selection import ColumnSelector

from sklearn.pipeline import make_pipeline

from sklearn.linear_model import LogisticRegression

iris = load_iris()

X = iris.data

y = iris.target

pipe1 = make_pipeline(ColumnSelector(cols=(0, 2)),

LogisticRegression())

pipe2 = make_pipeline(ColumnSelector(cols=(1, 2, 3)),

LogisticRegression())

sclf = StackingClassifier(classifiers=[pipe1, pipe2],

meta_classifier=LogisticRegression())

sclf.fit(X, y)

参考文档：
https://rasbt.github.io/mlxtend/user_guide/classifier/StackingClassifier/

18.2投票分类器（VotingClassifier）

投票分类器的原理是结合了多个不同的机器学习分类器，使用多数票或者平均预测概率（软票），预测类标签。这类分类器对一组相同表现的模型十分有用，同时可以平衡各自的弱点。投票分类又可进一步分为多数投票分类（Majority Class Labels）、加权平均概率（soft vote，软投票）。

18.2.1多数投票分类（MajorityVote Class）

多数投票分类的分类原则为预测标签不同时，按最多种类为最终分类；如果预测标签相同时，则按顺序，选择排在第1的标签为最终分类。举例如下：
预测类型的标签为该组学习器中相同最多的种类：例如给出的分类如下
分类器1 -> 标签1
分类器2 -> 标签1
分类器3 -> 标签2
投票分类器（voting=‘hard’）则该预测结果为‘标签1’。
在各个都只有一个的情况下，则按照顺序来，如下：
分类器1 -> 标签2
分类器2 -> 标签1
最终分类结果为“标签2”

18.2.1.1Iris数据集概述

首先，我们取得数据，下面这个链接中有数据的详细介绍，并可以下载数据集。https://archive.ics.uci.edu/ml/datasets/Iris
从数据的说明上，我们可以看到Iris有4个特征，3个类别。但是，我们为了数据的可视化，我们只保留2个特征（sepal length和petal length）。数据可视化代码如下：

%matplotlib inline
import pandas as pd
import matplotlib.pylab as plt
import numpy as np

# 加载Iris数据集作为DataFrame对象
df = pd.read_csv('http://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data', header=None)
 X = df.iloc[:, [0, 2]].values # 取出2个特征，并把它们用Numpy数组表示
 
plt.scatter(X[:50, 0], X[:50, 1],color='red', marker='o', label='setosa') # 前50个样本的散点图
plt.scatter(X[50:100, 0], X[50:100, 1],color='blue', marker='x', label='versicolor') # 中间50个样本的散点图
plt.scatter(X[100:, 0], X[100:, 1],color='green', marker='+', label='Virginica') # 后50个样本的散点图
plt.xlabel('petal length')
plt.ylabel('sepal length')
plt.legend(loc=2) # 把说明放在左上角，具体请参考官方文档
plt.show()

%matplotlib inline

import pandas as pd

import matplotlib.pylab as plt

import numpy as np

# 加载Iris数据集作为DataFrame对象

df = pd.read_csv('http://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data', header=None)

X = df.iloc[:, [0, 2]].values # 取出2个特征，并把它们用Numpy数组表示

plt.scatter(X[:50, 0], X[:50, 1],color='red', marker='o', label='setosa') # 前50个样本的散点图

plt.scatter(X[50:100, 0], X[50:100, 1],color='blue', marker='x', label='versicolor') # 中间50个样本的散点图

plt.scatter(X[100:, 0], X[100:, 1],color='green', marker='+', label='Virginica') # 后50个样本的散点图

plt.xlabel('petal length')

plt.ylabel('sepal length')

plt.legend(loc=2) # 把说明放在左上角，具体请参考官方文档

plt.show()

示例代码如下：

from sklearn import datasets
from sklearn import cross_validation
from sklearn.linear_model import LogisticRegression
from sklearn.naive_bayes import GaussianNB
from sklearn.ensemble import RandomForestClassifier
from sklearn.ensemble import VotingClassifier

iris = datasets.load_iris()
X, y = iris.data[:, 1:3], iris.target

clf1 = LogisticRegression(random_state=1)
clf2 = RandomForestClassifier(random_state=1)
clf3 = GaussianNB()

eclf = VotingClassifier(estimators=[('lr', clf1), ('rf', clf2), ('gnb', clf3)], voting='hard', weights=[2,1,2])

for clf, label in zip([clf1, clf2, clf3, eclf], ['Logistic Regression', 'Random Forest', 'naive Bayes', 'Ensemble']):
    scores = cross_validation.cross_val_score(clf, X, y, cv=5, scoring='accuracy')
    print("Accuracy: %0.2f (+/- %0.2f) [%s]" % (scores.mean(), scores.std(), label))

from sklearn import datasets

from sklearn import cross_validation

from sklearn.linear_model import LogisticRegression

from sklearn.naive_bayes import GaussianNB

from sklearn.ensemble import RandomForestClassifier

from sklearn.ensemble import VotingClassifier

iris = datasets.load_iris()

X, y = iris.data[:, 1:3], iris.target

clf1 = LogisticRegression(random_state=1)

clf2 = RandomForestClassifier(random_state=1)

clf3 = GaussianNB()

eclf = VotingClassifier(estimators=[('lr', clf1), ('rf', clf2), ('gnb', clf3)], voting='hard', weights=[2,1,2])

for clf, label in zip([clf1, clf2, clf3, eclf], ['Logistic Regression', 'Random Forest', 'naive Bayes', 'Ensemble']):

scores = cross_validation.cross_val_score(clf, X, y, cv=5, scoring='accuracy')

print("Accuracy: %0.2f (+/- %0.2f) [%s]" % (scores.mean(), scores.std(), label))

运行结果如下：
Accuracy: 0.90 (+/- 0.05) [Logistic Regression]
Accuracy: 0.93 (+/- 0.05) [Random Forest]
Accuracy: 0.91 (+/- 0.04) [naive Bayes]
Accuracy: 0.95 (+/- 0.05) [Ensemble]

18.2.2多数投票分类（MajorityVote Class）

相对于多数投票（hard voting），软投票返回预测概率值的总和最大的标签。可通过参数weights指定每个分类器的权重；若权重提供了，在计算时则会按照权重计算，然后取平均；标签则为概率最高的标签。
举例说明，假设有3个分类器，3个类，每个分类器的权重为：w1=1，w2=1，w3=1。如下表：

下面例子为线性SVM，决策树，K邻近分类器：

from sklearn import datasets  
from sklearn.tree import DecisionTreeClassifier  
from sklearn.neighbors import KNeighborsClassifier  
from sklearn.svm import SVC  
from itertools import product  
from sklearn.ensemble import VotingClassifier  

#Loading some example data  
iris = datasets.load_iris()  
X = iris.data[:, [0,2]]  
y = iris.target  

#Training classifiers  
clf1 = DecisionTreeClassifier(max_depth=4)  
clf2 = KNeighborsClassifier(n_neighbors=7)  
clf3 = SVC(kernel='rbf', probability=True)  
eclf = VotingClassifier(estimators=[('dt', clf1), ('knn', clf2), ('svc', clf3)], voting='soft', weights=[2,1,2])  

clf1 = clf1.fit(X,y)  
clf2 = clf2.fit(X,y)  
clf3 = clf3.fit(X,y)  
eclf = eclf.fit(X,y)


##这些分类器分类结果
x_min,x_max = X[:,0].min()-1,X[:,0].max()+1  
y_min,y_max = X[:,1].min()-1,X[:,1].max()+1  
xx,yy = np.meshgrid(np.arange(x_min,x_max,0.1),  
                    np.arange(y_min,y_max,0.1))  
f, axarr = plt.subplots(2, 2, sharex='col', sharey='row', figsize=(10, 8))  
for idx, clf, tt in zip(product([0, 1], [0, 1]),  
                        [clf1, clf2, clf3, eclf],  
                        ['Decision Tree (depth=4)', 'KNN (k=7)',  
                         'Kernel SVM', 'Soft Voting']):  
  
    Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])  
    Z = Z.reshape(xx.shape)  
  
    axarr[idx[0], idx[1]].contourf(xx, yy, Z, alpha=0.4)  
    axarr[idx[0], idx[1]].scatter(X[:, 0], X[:, 1], c=y, alpha=0.8)  
    axarr[idx[0], idx[1]].set_title(tt)  
plt.show()

from sklearn import datasets

from sklearn.tree import DecisionTreeClassifier

from sklearn.neighbors import KNeighborsClassifier

from sklearn.svm import SVC

from itertools import product

from sklearn.ensemble import VotingClassifier

#Loading some example data

iris = datasets.load_iris()

X = iris.data[:, [0,2]]

y = iris.target

#Training classifiers

clf1 = DecisionTreeClassifier(max_depth=4)

clf2 = KNeighborsClassifier(n_neighbors=7)

clf3 = SVC(kernel='rbf', probability=True)

eclf = VotingClassifier(estimators=[('dt', clf1), ('knn', clf2), ('svc', clf3)], voting='soft', weights=[2,1,2])

clf1 = clf1.fit(X,y)

clf2 = clf2.fit(X,y)

clf3 = clf3.fit(X,y)

eclf = eclf.fit(X,y)

##这些分类器分类结果

x_min,x_max = X[:,0].min()-1,X[:,0].max()+1

y_min,y_max = X[:,1].min()-1,X[:,1].max()+1

xx,yy = np.meshgrid(np.arange(x_min,x_max,0.1),

np.arange(y_min,y_max,0.1))

f, axarr = plt.subplots(2, 2, sharex='col', sharey='row', figsize=(10, 8))

for idx, clf, tt in zip(product([0, 1], [0, 1]),

[clf1, clf2, clf3, eclf],

['Decision Tree (depth=4)', 'KNN (k=7)',

'Kernel SVM', 'Soft Voting']):

Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])

Z = Z.reshape(xx.shape)

axarr[idx[0], idx[1]].contourf(xx, yy, Z, alpha=0.4)

axarr[idx[0], idx[1]].scatter(X[:, 0], X[:, 1], c=y, alpha=0.8)

axarr[idx[0], idx[1]].set_title(tt)

plt.show()

18.3自适应分类器（Adaboost）

Adaboost是一种迭代算法，其核心思想是针对同一个训练集训练不同的分类器(弱分类器)，然后把这些弱分类器集合起来，构成一个更强的最终分类器(强分类器)。其算法本身是通过改变数据分布来实现的，它根据每次训练集之中每个样本的分类是否正确，以及上次的总体分类的准确率，来确定每个样本的权值。将修改过权值的新数据集送给下层分类器进行训练，最后将每次训练得到的分类器最后融合起来，作为最后的决策分类器。使用adaboost分类器可以排除一些不必要的训练数据特征，并放在关键的训练数据上面。
下面的例子展示了AdaBoost算法拟合100个弱学习器

from sklearn.model_selection import cross_val_score  
from sklearn.datasets import load_iris  
from sklearn.ensemble import AdaBoostClassifier  

iris = load_iris()  
clf = AdaBoostClassifier(n_estimators=100)  
scores = cross_val_score(clf, iris.data, iris.target)  
scores.mean()

from sklearn.model_selection import cross_val_score

from sklearn.datasets import load_iris

from sklearn.ensemble import AdaBoostClassifier

iris = load_iris()

clf = AdaBoostClassifier(n_estimators=100)

scores = cross_val_score(clf, iris.data, iris.target)

scores.mean()

输出结果为：
0.95996732026143794

18.4 Xgboost简介

18.4.1简介

Xgboost是很多CART回归树集成，CART树以基尼系数为划分依据。回归树的样本输出是数值的形式，比如给某人发放房屋贷款的数额就是具体的数值，可以是0到120万元之间的任意值。那么，这时候你就没法用上述的信息增益、信息增益率、基尼系数来判定树的节点分裂了，你就会采用新的方式，预测误差，常用的有均方误差、对数误差等。而且节点不再是类别，是数值（预测值），那么怎么确定呢，有的是节点内样本均值，有的是最优化算出来的比如Xgboost。
xgboot特点：
1、w是最优化求出来的
2、使用正则化防止过拟合的技术
3、支持分布式、并行化，树之间没有很强的上下依赖关系
4、支持GPU
下图就是一个CART的例子，CART会把输入根据属性分配到各个也子节点上，而每个叶子节点上面会对应一个分数值。下面的例子是预测一个人是否喜欢电脑游戏。将叶子节点表示为分数之后，可以做很多事情，比如概率预测，排序等等。

一个CART往往过于简单，而无法有效的进行预测，因此更加高效的是使用多个CART进行融合，使用集成的方法提升预测效率：

假设有两颗回归树，则两棵树融合后的预测结果如上图。
xgboost涉及参数较多，具体使用可参考：https://cloud.tencent.com/developer/article/1111048
https://blog.csdn.net/han_xiaoyang/article/details/52665396

18.4.2 xgboost 实例

主要目的：利用多种预测方法，对房价进行预测
数据结构：
数据探索与预处理：
创建及优化模型：

18.4.2.1数据探索及数据预处理

（1）导入需要的库，并导入数据，查看前五行样本数据，参考文档：

#import some necessary librairies

import numpy as np # linear algebra
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)
%matplotlib inline
import matplotlib.pyplot as plt  # Matlab-style plotting
import seaborn as sns
color = sns.color_palette()
sns.set_style('darkgrid')
import warnings
def ignore_warn(*args, **kwargs):
    pass
warnings.warn = ignore_warn #ignore annoying warning (from sklearn and seaborn)


from scipy import stats
from scipy.stats import norm, skew #for some statistics


pd.set_option('display.float_format', lambda x: '{:.3f}'.format(x)) #Limiting floats output to 3 decimal points


from subprocess import check_output
print(check_output(["ls", "../data"]).decode("utf8")) #check the files available in the directory


#Now let's import and put the train and test datasets in  pandas dataframe

train = pd.read_csv('../data/house_train.csv')
test = pd.read_csv('../data/house_test.csv')
##display the first five rows of the train dataset.
train.head(5)
 
#check the numbers of samples and features
print("The train data size before dropping Id feature is : {} ".format(train.shape))
print("The test data size before dropping Id feature is : {} ".format(test.shape))

#Save the 'Id' column
train_ID = train['Id']
test_ID = test['Id']

#对原数组直接删除ID列，inplace = True为直接删除原数组.
train.drop("Id", axis = 1, inplace = True)
test.drop("Id", axis = 1, inplace = True)

#check again the data size after dropping the 'Id' variable
print("\nThe train data size after dropping Id feature is : {} ".format(train.shape)) 
print("The test data size after dropping Id feature is : {} ".format(test.shape))

#import some necessary librairies

import numpy as np # linear algebra

import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)

%matplotlib inline

import matplotlib.pyplot as plt # Matlab-style plotting

import seaborn as sns

color = sns.color_palette()

sns.set_style('darkgrid')

import warnings

def ignore_warn(*args, **kwargs):

pass

warnings.warn = ignore_warn #ignore annoying warning (from sklearn and seaborn)

from scipy import stats

from scipy.stats import norm, skew #for some statistics

pd.set_option('display.float_format', lambda x: '{:.3f}'.format(x)) #Limiting floats output to 3 decimal points

from subprocess import check_output

print(check_output(["ls", "../data"]).decode("utf8")) #check the files available in the directory

#Now let's import and put the train and test datasets in pandas dataframe

train = pd.read_csv('../data/house_train.csv')

test = pd.read_csv('../data/house_test.csv')

##display the first five rows of the train dataset.

train.head(5)

#check the numbers of samples and features

print("The train data size before dropping Id feature is : {} ".format(train.shape))

print("The test data size before dropping Id feature is : {} ".format(test.shape))

#Save the 'Id' column

train_ID = train['Id']

test_ID = test['Id']

#对原数组直接删除ID列，inplace = True为直接删除原数组.

train.drop("Id", axis = 1, inplace = True)

test.drop("Id", axis = 1, inplace = True)

#check again the data size after dropping the 'Id' variable

print("\nThe train data size after dropping Id feature is : {} ".format(train.shape))

print("The test data size after dropping Id feature is : {} ".format(test.shape))

The train data size before dropping Id feature is : (1460, 81)
The test data size before dropping Id feature is : (1459, 80)

The train data size after dropping Id feature is : (1460, 80)
The test data size after dropping Id feature is : (1459, 79)

（2）探索孤立点

fig, ax = plt.subplots()
ax.scatter(x = train['GrLivArea'], y = train['SalePrice'])
plt.ylabel('SalePrice', fontsize=13)
plt.xlabel('GrLivArea', fontsize=13)
plt.show()

fig, ax = plt.subplots()

ax.scatter(x = train['GrLivArea'], y = train['SalePrice'])

plt.ylabel('SalePrice', fontsize=13)

plt.xlabel('GrLivArea', fontsize=13)

plt.show()

（3）删除一些孤立点
删除房屋销售价（SalePr ice）小于300000（美元）并且居住面积（GrLivArea）平方英尺大于40003的记录。主要是上图中右下边这几个点。

#Deleting outliers
train = train.drop(train[(train['GrLivArea']>4000) & (train['SalePrice']<300000)].index)

#Check the graphic again
fig, ax = plt.subplots()
ax.scatter(train['GrLivArea'], train['SalePrice'])
plt.ylabel('SalePrice', fontsize=13)
plt.xlabel('GrLivArea', fontsize=13)
plt.show()

#Deleting outliers

train = train.drop(train[(train['GrLivArea']>4000) & (train['SalePrice']<300000)].index)

#Check the graphic again

fig, ax = plt.subplots()

ax.scatter(train['GrLivArea'], train['SalePrice'])

plt.ylabel('SalePrice', fontsize=13)

plt.xlabel('GrLivArea', fontsize=13)

plt.show()

（4）探索房价分布情况
分析目标变量（房价），画出房价分布图及QQ图，QQ图就是分位数图示法（Quantile Quantile Plot，Q-Q图主要用于检验数据分布的相似性，如果要利用Q-Q图来对数据进行正态分布的检验，则可以令x轴为正态分布的分位数，y轴为样本分位数，如果这两者构成的点分布在一条直线上，就证明样本数据与正态分布存在线性相关性，即服从正态分布。

sns.distplot(train['SalePrice'] , fit=norm);

# Get the fitted parameters used by the function
(mu, sigma) = norm.fit(train['SalePrice'])
print( '\n mu = {:.2f} and sigma = {:.2f}\n'.format(mu, sigma))

#Now plot the distribution
plt.legend(['Normal dist. ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)],
            loc='best')
plt.ylabel('Frequency')
plt.title('SalePrice distribution')

#Get also the QQ-plot
fig = plt.figure()
res = stats.probplot(train['SalePrice'], plot=plt)
plt.show()

sns.distplot(train['SalePrice'] , fit=norm);

# Get the fitted parameters used by the function

(mu, sigma) = norm.fit(train['SalePrice'])

print( '\n mu = {:.2f} and sigma = {:.2f}\n'.format(mu, sigma))

#Now plot the distribution

plt.legend(['Normal dist. ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)],

loc='best')

plt.ylabel('Frequency')

plt.title('SalePrice distribution')

#Get also the QQ-plot

fig = plt.figure()

res = stats.probplot(train['SalePrice'], plot=plt)

plt.show()

由上图可知，目标变量是右倾斜的。由于（线性）模型喜欢正态分布的数据，我们需要对房价特征进行转换，使其接近正态分布。

（5）对房价特征进行log转换

#We use the numpy fuction log1p which  applies log(1+x) to all elements of the column
train["SalePrice"] = np.log1p(train["SalePrice"])

#Check the new distribution 
sns.distplot(train['SalePrice'] , fit=norm);

# Get the fitted parameters used by the function
(mu, sigma) = norm.fit(train['SalePrice'])
print( '\n mu = {:.2f} and sigma = {:.2f}\n'.format(mu, sigma))

#Now plot the distribution
plt.legend(['Normal dist. ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)],
            loc='best')
plt.ylabel('Frequency')
plt.title('SalePrice distribution')
#Get also the QQ-plot
fig = plt.figure()
res = stats.probplot(train['SalePrice'], plot=plt)
plt.show()

#We use the numpy fuction log1p which applies log(1+x) to all elements of the column

train["SalePrice"] = np.log1p(train["SalePrice"])

#Check the new distribution

sns.distplot(train['SalePrice'] , fit=norm);

# Get the fitted parameters used by the function

(mu, sigma) = norm.fit(train['SalePrice'])

print( '\n mu = {:.2f} and sigma = {:.2f}\n'.format(mu, sigma))

#Now plot the distribution

plt.legend(['Normal dist. ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)],

loc='best')

plt.ylabel('Frequency')

plt.title('SalePrice distribution')

#Get also the QQ-plot

fig = plt.figure()

res = stats.probplot(train['SalePrice'], plot=plt)

plt.show()

现在纠正了偏差，数据看起来更正常分布。
（6）连接训练数据和测试数据
为便于统一处理，我们需要把训练数据、测试数据集成在一起。对房价数据不做处理。

ntrain = train.shape[0]
ntest = test.shape[0]
y_train = train.SalePrice.values
#为使合并后的索引正常，加上reset_index
all_data = pd.concat((train, test)).reset_index(drop=True)
all_data.drop(['SalePrice'], axis=1, inplace=True)
print("all_data size is : {}".format(all_data.shape))

ntrain = train.shape[0]

ntest = test.shape[0]

y_train = train.SalePrice.values

#为使合并后的索引正常，加上reset_index

all_data = pd.concat((train, test)).reset_index(drop=True)

all_data.drop(['SalePrice'], axis=1, inplace=True)

print("all_data size is : {}".format(all_data.shape))

运行结果：
all_data size is : (2917, 79)
（7）查看缺失数据情况
以下我们查看各特征的缺失率

all_data_na = (all_data.isnull().sum() / len(all_data)) * 100
all_data_na = all_data_na.drop(all_data_na[all_data_na == 0].index).sort_values(ascending=False)[:30]
missing_data = pd.DataFrame({'Missing Ratio' :all_data_na})
missing_data.head(20)

all_data_na = (all_data.isnull().sum() / len(all_data)) * 100

all_data_na = all_data_na.drop(all_data_na[all_data_na == 0].index).sort_values(ascending=False)[:30]

missing_data = pd.DataFrame({'Missing Ratio' :all_data_na})

missing_data.head(20)

运行结果
Missing Ratio
PoolQC 99.691
MiscFeature 96.4
Alley 93.212
Fence 80.425
FireplaceQu 48.68
LotFrontage 16.661
GarageQual 5.451
GarageCond 5.451
GarageFinish 5.451
GarageYrBlt 5.451
GarageType 5.382
BsmtExposure 2.811
BsmtCond 2.811
BsmtQual 2.777
BsmtFinType2 2.743
BsmtFinType1 2.708
MasVnrType 0.823
MasVnrArea 0.788
MSZoning 0.137
BsmtFullBath 0.069
可视化这些数据

f, ax = plt.subplots(figsize=(15, 12))
plt.xticks(rotation='90')
sns.barplot(x=all_data_na.index, y=all_data_na)
plt.xlabel('Features', fontsize=15)
plt.ylabel('Percent of missing values', fontsize=15)
plt.title('Percent missing data by feature', fontsize=15)

f, ax = plt.subplots(figsize=(15, 12))

plt.xticks(rotation='90')

sns.barplot(x=all_data_na.index, y=all_data_na)

plt.xlabel('Features', fontsize=15)

plt.ylabel('Percent of missing values', fontsize=15)

plt.title('Percent missing data by feature', fontsize=15)

（8）查看数据的相关性

#Correlation map to see how features are correlated with SalePrice
corrmat = train.corr()
plt.subplots(figsize=(12,9))
sns.heatmap(corrmat, vmax=0.9, square=True)

#Correlation map to see how features are correlated with SalePrice

corrmat = train.corr()

plt.subplots(figsize=(12,9))

sns.heatmap(corrmat, vmax=0.9, square=True)

颜色越深，表示相关性越强。
（9）填充缺失值
以下我们对存在缺失值的特征分别进行处理，
PoolQC：数据描述表示NA表示“无池”。这是有道理的，因为缺失值的比例很大（+ 99％），而且大多数房屋一般都没有游泳池。这里我们把缺失值改为None。

all_data["PoolQC"] = all_data["PoolQC"].fillna("None")

1	all_data["PoolQC"] = all_data["PoolQC"].fillna("None")

MiscFeature：数据描述表示NA表示“没有杂项功能，这里把缺失值改为None
all_data["MiscFeature"] = all_data["MiscFeature"].fillna("None")
类似把特征Alley 、Fence、FireplaceQu、GarageType, GarageFinish, GarageQual and GarageCond,进行相同处理，这些特征的缺失率都比较高。

for col in ('Alley','Fence','FireplaceQu','GarageType', 'GarageFinish', 'GarageQual', 'GarageCond'):
    all_data[col] = all_data[col].fillna('None')

1 2	for col in ('Alley','Fence','FireplaceQu','GarageType', 'GarageFinish', 'GarageQual', 'GarageCond'): all_data[col] = all_data[col].fillna('None')

LotFrontage（与街道连接的线性脚）：由于连接到房产的每条街道的区域很可能与其附近的其他房屋有相似的区域，我们可以通过邻域的中位数LotFrontage填写缺失值。

#Group by neighborhood and fill in missing value by the median LotFrontage of all the neighborhood
all_data["LotFrontage"] = all_data.groupby("Neighborhood")["LotFrontage"].transform(
    lambda x: x.fillna(x.median()))

#Group by neighborhood and fill in missing value by the median LotFrontage of all the neighborhood

all_data["LotFrontage"] = all_data.groupby("Neighborhood")["LotFrontage"].transform(

lambda x: x.fillna(x.median()))

GarageYrBlt，GarageArea和GarageCars：用0代替缺失数据（因为没有车库=这样的车库没有车辆。）

for col in ('GarageYrBlt', 'GarageArea', 'GarageCars'):
    all_data[col] = all_data[col].fillna(0)

1 2	for col in ('GarageYrBlt', 'GarageArea', 'GarageCars'): all_data[col] = all_data[col].fillna(0)

对BsmtFinSF1, BsmtFinSF2, BsmtUnfSF, TotalBsmtSF, BsmtFullBath and BsmtHalfBath, MasVnrArea做同样处理。

for col in ('BsmtFinSF1', 'BsmtFinSF2', 'BsmtUnfSF','TotalBsmtSF', 'BsmtFullBath', 'BsmtHalfBath','MasVnrArea'):
    all_data[col] = all_data[col].fillna(0)

1 2	for col in ('BsmtFinSF1', 'BsmtFinSF2', 'BsmtUnfSF','TotalBsmtSF', 'BsmtFullBath', 'BsmtHalfBath','MasVnrArea'): all_data[col] = all_data[col].fillna(0)

对BsmtQual', 'BsmtCond', 'BsmtExposure', 'BsmtFinType1', 'BsmtFinType2','MasVnrType特征的缺失值改为None。

for col in ('BsmtQual', 'BsmtCond', 'BsmtExposure', 'BsmtFinType1', 'BsmtFinType2','MasVnrType'):
    all_data[col] = all_data[col].fillna('None')

1 2	for col in ('BsmtQual', 'BsmtCond', 'BsmtExposure', 'BsmtFinType1', 'BsmtFinType2','MasVnrType'): all_data[col] = all_data[col].fillna('None')

MSZoning（一般分区分类）：'RL'是迄今为止最常见的值。所以我们可以用'RL'来填补缺失值,取频度最大的那个数据，采用mode()[0]的格式。

all_data['MSZoning'] = all_data['MSZoning'].fillna(all_data['MSZoning'].mode()[0])

1	all_data['MSZoning'] = all_data['MSZoning'].fillna(all_data['MSZoning'].mode()[0])

对特征Electrical,KitchenQual,Exterior1st,Exterior2nd,SaleType做相同处理。

for col in ('Electrical','KitchenQual','Exterior1st','Exterior2nd','SaleType'):
    all_data[col] = all_data[col].fillna(all_data[col].mode()[0])

1 2	for col in ('Electrical','KitchenQual','Exterior1st','Exterior2nd','SaleType'): all_data[col] = all_data[col].fillna(all_data[col].mode()[0])

Utilities：对于此分类功能，所有记录都是“AllPub”，除了一个“NoSeWa”和2个NA。由于带有“NoSewa”的房子位于训练集中，因此该功能无助于预测建模。然后我们可以安全地删除它。

all_data = all_data.drop(['Utilities'], axis=1)
Functional：数据描述说NA意味着典型
all_data["Functional"] = all_data["Functional"].fillna("Typ")

MSSubClass：Na很可能意味着没有建筑类。 我们可以用None替换缺失值

all_data = all_data.drop(['Utilities'], axis=1)

Functional：数据描述说NA意味着典型

all_data["Functional"] = all_data["Functional"].fillna("Typ")

MSSubClass：Na很可能意味着没有建筑类。我们可以用None替换缺失值

all_data['MSSubClass'] = all_data['MSSubClass'].fillna("None")

1	all_data['MSSubClass'] = all_data['MSSubClass'].fillna("None")

通过以上缺失值的填充，我们再看一下缺失情况

#Check remaining missing values if any 
all_data_na = (all_data.isnull().sum() / len(all_data)) * 100
all_data_na = all_data_na.drop(all_data_na[all_data_na == 0].index).sort_values(ascending=False)
missing_data = pd.DataFrame({'Missing Ratio' :all_data_na})
missing_data.head()

#Check remaining missing values if any

all_data_na = (all_data.isnull().sum() / len(all_data)) * 100

all_data_na = all_data_na.drop(all_data_na[all_data_na == 0].index).sort_values(ascending=False)

missing_data = pd.DataFrame({'Missing Ratio' :all_data_na})

missing_data.head()

运行结果为
没有缺失值特征了！
（10）对一些分类的数值特征转换。
这里使用LabelEncoder 对不连续的数字或者文本进行编号。

#MSSubClass=The building class
all_data['MSSubClass'] = all_data['MSSubClass'].apply(str)


#Changing OverallCond into a categorical variable
all_data['OverallCond'] = all_data['OverallCond'].astype(str)


#Year and month sold are transformed into categorical features.
all_data['YrSold'] = all_data['YrSold'].astype(str)
all_data['MoSold'] = all_data['MoSold'].astype(str)

#MSSubClass=The building class

all_data['MSSubClass'] = all_data['MSSubClass'].apply(str)

#Changing OverallCond into a categorical variable

all_data['OverallCond'] = all_data['OverallCond'].astype(str)

#Year and month sold are transformed into categorical features.

all_data['YrSold'] = all_data['YrSold'].astype(str)

all_data['MoSold'] = all_data['MoSold'].astype(str)

【说明】
apply(str)与astype(str)功能相同，一般前者快于后者。

（11）对类别特征进行编码
这里使用LabelEncoder 对不连续的数字或者文本进行编号。

from sklearn.preprocessing import LabelEncoder
cols = ('FireplaceQu', 'BsmtQual', 'BsmtCond', 'GarageQual', 'GarageCond', 
        'ExterQual', 'ExterCond','HeatingQC', 'PoolQC', 'KitchenQual', 'BsmtFinType1', 
        'BsmtFinType2', 'Functional', 'Fence', 'BsmtExposure', 'GarageFinish', 'LandSlope',
        'LotShape', 'PavedDrive', 'Street', 'Alley', 'CentralAir', 'MSSubClass', 'OverallCond', 
        'YrSold', 'MoSold')
# process columns, apply LabelEncoder to categorical features
for c in cols:
    lbl = LabelEncoder() 
    lbl.fit(list(all_data[c].values)) 
    all_data[c] = lbl.transform(list(all_data[c].values))

# shape        
print('Shape all_data: {}'.format(all_data.shape))

from sklearn.preprocessing import LabelEncoder

cols = ('FireplaceQu', 'BsmtQual', 'BsmtCond', 'GarageQual', 'GarageCond',

'ExterQual', 'ExterCond','HeatingQC', 'PoolQC', 'KitchenQual', 'BsmtFinType1',

'BsmtFinType2', 'Functional', 'Fence', 'BsmtExposure', 'GarageFinish', 'LandSlope',

'LotShape', 'PavedDrive', 'Street', 'Alley', 'CentralAir', 'MSSubClass', 'OverallCond',

'YrSold', 'MoSold')

# process columns, apply LabelEncoder to categorical features

for c in cols:

lbl = LabelEncoder()

lbl.fit(list(all_data[c].values))

all_data[c] = lbl.transform(list(all_data[c].values))

# shape

print('Shape all_data: {}'.format(all_data.shape))

（12）增加特征
由于区域相关特征对于确定房价非常重要，这里将增加了一个特征，即每个房屋的地下室总面积，一楼和二楼面积之和

#Adding total sqfootage feature 
all_data['TotalSF'] = all_data['TotalBsmtSF'] + all_data['1stFlrSF'] + all_data['2ndFlrSF']

1 2	#Adding total sqfootage feature all_data['TotalSF'] = all_data['TotalBsmtSF'] + all_data['1stFlrSF'] + all_data['2ndFlrSF']

（13）查看一些特征的歪斜程度
这里主要对一些数值型特征进行分析

numeric_feats = all_data.dtypes[all_data.dtypes != "object"].index

# Check the skew of all numerical features
skewed_feats = all_data[numeric_feats].apply(lambda x: skew(x.dropna())).sort_values(ascending=False)
print("\nSkew in numerical features: \n")
skewness = pd.DataFrame({'Skew' :skewed_feats})
skewness.head(10)

numeric_feats = all_data.dtypes[all_data.dtypes != "object"].index

# Check the skew of all numerical features

skewed_feats = all_data[numeric_feats].apply(lambda x: skew(x.dropna())).sort_values(ascending=False)

print("\nSkew in numerical features: \n")

skewness = pd.DataFrame({'Skew' :skewed_feats})

skewness.head(10)

运行结果

（14）对这些倾斜特征进行box-cox变换
Box-Cox变换是统计建模中常用的一种数据变换，用于连续的响应变量不满足正态分布的情况。Box-Cox变换使线性回归模型满足线性性、独立性、方差齐性以及正态性的同时，又不丢失信息。
使用Box-Cox变换族一般都可以保证将数据进行成功的正态变换，但在二分变量或较少水平的等级变量的情况下，不能成功进行转换，此时，我们可以考虑使用广义线性模型，如LOGUSTICS模型、Johnson转换等.
Box-Cox变换后，残差可以更好的满足正态性、独立性等假设前提，降低了伪回归的概率

其中确定λ是关键，如何确定λ？一般采用最大似然估计来求。该值一般在[-5,5]之间。
一般y分布为左偏(左边比较陡峭)，则λ=0；如果y为右偏，取λ>0。

skewness = skewness[abs(skewness) > 0.75]
print("There are {} skewed numerical features to Box Cox transform".format(skewness.shape[0]))

from scipy.special import boxcox1p
skewed_features = skewness.index
lam = 0.15
for feat in skewed_features:
    #all_data[feat] += 1
    all_data[feat] = boxcox1p(all_data[feat], lam)

skewness = skewness[abs(skewness) > 0.75]

print("There are {} skewed numerical features to Box Cox transform".format(skewness.shape[0]))

from scipy.special import boxcox1p

skewed_features = skewness.index

lam = 0.15

for feat in skewed_features:

#all_data[feat] += 1

all_data[feat] = boxcox1p(all_data[feat], lam)

（15）得到虚拟分类特征，对标称类别转换为one-hot编码。

all_data = pd.get_dummies(all_data)
print(all_data.shape)

1 2	all_data = pd.get_dummies(all_data) print(all_data.shape)

运行结果：
(2917, 220)
（16）得到新的训练集、测试集

train = all_data[:ntrain]
test = all_data[ntrain:]

1 2	train = all_data[:ntrain] test = all_data[ntrain:]

查看他们的维度

print(train.shape)
print(test.shape)

1 2	print(train.shape) print(test.shape)

运行结果
(1458, 220)
(1459, 220)
至此，数据探索及预处理就基本完成，接下来开始创建模型。

18.4.2.2创建模型

（1）导入需要的库

from sklearn.linear_model import ElasticNet, Lasso,  BayesianRidge, LassoLarsIC
from sklearn.ensemble import RandomForestRegressor,  GradientBoostingRegressor
from sklearn.kernel_ridge import KernelRidge
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import RobustScaler
from sklearn.base import BaseEstimator, TransformerMixin, RegressorMixin, clone
from sklearn.model_selection import KFold, cross_val_score, train_test_split
from sklearn.metrics import mean_squared_error
import xgboost as xgb
import lightgbm as lgb

from sklearn.linear_model import ElasticNet, Lasso, BayesianRidge, LassoLarsIC

from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor

from sklearn.kernel_ridge import KernelRidge

from sklearn.pipeline import make_pipeline

from sklearn.preprocessing import RobustScaler

from sklearn.base import BaseEstimator, TransformerMixin, RegressorMixin, clone

from sklearn.model_selection import KFold, cross_val_score, train_test_split

from sklearn.metrics import mean_squared_error

import xgboost as xgb

import lightgbm as lgb

（2）使用K折交叉验证，其中k=5
K折交叉验证简介：
将数据集平均分割成K个等份
使用1份数据作为测试数据，其余作为训练数据
计算测试准确率
使用不同的测试集，重复2、3步骤
对测试准确率做平均，作为对未知数据预测准确率的估计
如下图

#Validation function
n_folds = 5

def rmsle_cv(model):
    kf = KFold(n_folds, shuffle=True, random_state=42).get_n_splits(train.values)
    rmse= np.sqrt(-cross_val_score(model, train.values, y_train, scoring="neg_mean_squared_error", cv = kf))
    return(rmse)

#Validation function

n_folds = 5

def rmsle_cv(model):

kf = KFold(n_folds, shuffle=True, random_state=42).get_n_splits(train.values)

rmse= np.sqrt(-cross_val_score(model, train.values, y_train, scoring="neg_mean_squared_error", cv = kf))

return(rmse)

【说明】
# 这里的cross_val_score将交叉验证的整个过程连接起来，不用再进行手动的分割数据
# cv参数用于规定将原始数据分成多少份
# scoring：该参数来控制它们对 estimators evaluated （评估的估计量）应用的指标。
对分类模型，该值可以为‘accuracy’ 或‘f1’；对回归模型，可以为‘explained_variance’， ‘neg_mean_squared_error’等。对性能来说，越小越好，最佳为0.
#多种评估指标：回归模型
均方差（MSE）:

均方误差对数（MSLE）：

平均绝对误差（MAE）:

更多模型评估指标，可参考：http://sklearn.apachecn.org/cn/0.19.0/modules/model_evaluation.html
（3）使用模型
在模型选择上，“没有免费的午餐”。为了比较模型间的性能，这里使用多种模型。首先使用3 种回归模型（Linear Regression，Lasso，Ridge。
使用LASSO Regression :
该模型可能对异常值非常敏感。所以我们需要让它们更加健壮。为此，我们在管道上使用sklearn的Robustscaler（）方法。

lasso = make_pipeline(RobustScaler(), Lasso(alpha =0.0005, random_state=1))

1	lasso = make_pipeline(RobustScaler(), Lasso(alpha =0.0005, random_state=1))

使用Elastic Net Regression :
增益对异常值有所增强，故这里也使用Robustscaler（）方法

ENet = make_pipeline(RobustScaler(), ElasticNet(alpha=0.0005, l1_ratio=.9, random_state=3))

1	ENet = make_pipeline(RobustScaler(), ElasticNet(alpha=0.0005, l1_ratio=.9, random_state=3))

【ElasticNet简介】
ElasticNet 是一种使用L1和L2先验作为正则化矩阵的线性回归模型.这种组合用于只有很少的权重非零的稀疏模型，比如:class:Lasso, 但是又能保持:class:Ridge 的正则化属性.我们可以使用ρ（ρ的具体位置，请看下面这个表达式）参数来调节L1和L2的凸组合(一类特殊的线性组合)。
当多个特征和另一个特征相关的时候弹性网络非常有用。Lasso 倾向于随机选择其中一个，而弹性网络更倾向于选择两个.
在实践中，Lasso 和 Ridge 之间权衡的一个优势是它允许在循环过程（Under rotate）中继承 Ridge 的稳定性.
弹性网络的目标函数是最小化：

ElasticNetCV 可以通过交叉验证来用来设置参数:
alpha (α)，l1_ratio (ρ)

使用Kernel Ridge Regression :

KRR = KernelRidge(alpha=0.6, kernel='polynomial', degree=2, coef0=2.5)

1	KRR = KernelRidge(alpha=0.6, kernel='polynomial', degree=2, coef0=2.5)

使用Gradient Boosting Regression :
其对异常值具有鲁棒性

GBoost = GradientBoostingRegressor(n_estimators=3000, learning_rate=0.05,
                                   max_depth=4, max_features='sqrt',
                                   min_samples_leaf=15, min_samples_split=10, 
                                   loss='huber', random_state =5)

GBoost = GradientBoostingRegressor(n_estimators=3000, learning_rate=0.05,

max_depth=4, max_features='sqrt',

min_samples_leaf=15, min_samples_split=10,

loss='huber', random_state =5)

使用XGBoost :

model_xgb = xgb.XGBRegressor(colsample_bytree=0.4603, gamma=0.0468, 
                             learning_rate=0.05, max_depth=3, 
                             min_child_weight=1.7817, n_estimators=2200,
                             reg_alpha=0.4640, reg_lambda=0.8571,
                             subsample=0.5213, silent=1,
                             random_state =7, nthread = -1)

model_xgb = xgb.XGBRegressor(colsample_bytree=0.4603, gamma=0.0468,

learning_rate=0.05, max_depth=3,

min_child_weight=1.7817, n_estimators=2200,

reg_alpha=0.4640, reg_lambda=0.8571,

subsample=0.5213, silent=1,

random_state =7, nthread = -1)

使用LightGBM

model_lgb = lgb.LGBMRegressor(objective='regression',num_leaves=5,
                              learning_rate=0.05, n_estimators=720,
                              max_bin = 55, bagging_fraction = 0.8,
                              bagging_freq = 5, feature_fraction = 0.2319,
                              feature_fraction_seed=9, bagging_seed=9,
                              min_data_in_leaf =6, min_sum_hessian_in_leaf = 11)

model_lgb = lgb.LGBMRegressor(objective='regression',num_leaves=5,

learning_rate=0.05, n_estimators=720,

max_bin = 55, bagging_fraction = 0.8,

bagging_freq = 5, feature_fraction = 0.2319,

feature_fraction_seed=9, bagging_seed=9,

min_data_in_leaf =6, min_sum_hessian_in_leaf = 11)

18.4.2.3看各模型的性能

让我们通过评估交叉验证rmsle错误来了解这些基本模型如何对数据执行

score = rmsle_cv(lasso)
print("\nLasso score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))

score = rmsle_cv(ENet)
print("ElasticNet score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))

score = rmsle_cv(KRR)
print("Kernel Ridge score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))

score = rmsle_cv(GBoost)
print("Gradient Boosting score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))

score = rmsle_cv(model_xgb)
print("Xgboost score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))

score = rmsle_cv(model_lgb)
print("LGBM score: {:.4f} ({:.4f})\n" .format(score.mean(), score.std()))

score = rmsle_cv(lasso)

print("\nLasso score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))

score = rmsle_cv(ENet)

print("ElasticNet score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))

score = rmsle_cv(KRR)

print("Kernel Ridge score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))

score = rmsle_cv(GBoost)

print("Gradient Boosting score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))

score = rmsle_cv(model_xgb)

print("Xgboost score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))

score = rmsle_cv(model_lgb)

print("LGBM score: {:.4f} ({:.4f})\n" .format(score.mean(), score.std()))

运行结果
Lasso score: 0.0032 (0.0002)
ElasticNet score: 0.0031 (0.0002)
Kernel Ridge score: 0.0034 (0.0004)
Gradient Boosting score: 0.0026 (0.0002)
Xgboost score: 0.0086 (0.0003)
LGBM score: 0.0026 (0.0002)

18.4.2.4堆叠（或集成）各模型

我们从这种平均基本模型的简单方法开始。我们构建了一个新类来扩展scikit-learn与我们的模型，并且还包括封装和代码重用（继承）
（1）对各模型求平均值

class AveragingModels(BaseEstimator, RegressorMixin, TransformerMixin):
    def __init__(self, models):
        self.models = models
        
    # we define clones of the original models to fit the data in
    def fit(self, X, y):
        self.models_ = [clone(x) for x in self.models]
        
        # Train cloned base models
        for model in self.models_:
            model.fit(X, y)

        return self
    
    #Now we do the predictions for cloned models and average them
    def predict(self, X):
        predictions = np.column_stack([
            model.predict(X) for model in self.models_
        ])
        return np.mean(predictions, axis=1)

class AveragingModels(BaseEstimator, RegressorMixin, TransformerMixin):

def __init__(self, models):

self.models = models

# we define clones of the original models to fit the data in

def fit(self, X, y):

self.models_ = [clone(x) for x in self.models]

# Train cloned base models

for model in self.models_:

model.fit(X, y)

return self

#Now we do the predictions for cloned models and average them

def predict(self, X):

predictions = np.column_stack([

model.predict(X) for model in self.models_

])

return np.mean(predictions, axis=1)

（2）看集成后模型性能
这里集成ENet, GBoost, KRR, lasso这4种基本模型

averaged_models = AveragingModels(models = (ENet, GBoost, KRR, lasso))

score = rmsle_cv(averaged_models)
print(" Averaged base models score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))

averaged_models = AveragingModels(models = (ENet, GBoost, KRR, lasso))

score = rmsle_cv(averaged_models)

print(" Averaged base models score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))

运行结果
Averaged base models score: 0.0027 (0.0002)
不错！从初步结果来看，上面我们采用的最简单的堆叠方法确实提高了分数。这将鼓励我们进一步探索不那么简单的堆叠方法。

18.4.2.5添加元模型

这种集成方法的基本思想如下：
在这种方法中，我们在平均基础模型上添加元模型，并使用这些基础模型的折叠后预测来训练我们的元模型。
培训部分的程序可以描述如下：
1）将整个训练集分成两个不相交的集（这里是train和holdout）
2）在第一部分（即train）训练几个基础模型
3）在第二部分（即holdout）测试这些基础模型
4）使用来自3）的预测（称为折叠外预测）作为输入，并使用正确的响应（目标变量）作为输出来训练更高级别的元模型。
前三个步骤是迭代完成的。如果我们采用5-折（5-fold）堆叠，我们首先将训练数据分成5折（fold）。然后我们将进行5次迭代。在每次迭代中，我们训练每个基础模型4份并在剩余的一份数据上进行预测。
因此，经过5次迭代后，我们将确保使用整个数据进行折叠后预测，然后我们将使用这些预测作为新特征来训练第4步中的元模型。
对于预测部分，我们对测试数据上所有基础模型的预测进行平均，并将它们用作元特征，最终预测是使用元模型完成的。

以上步骤如下图所示：

数据生成过程的动态效果图：

有关stacking的介绍可参考：
http://blog.kaggle.com/2017/06/15/stacking-made-easy-an-introduction-to-stacknet-by-competitions-grandmaster-marios-michailidis-kazanova/
此外，下图从另一个角度来说明stacking的原理。
改为如下图（主要把上部分的列标题都为Model1）

如何理解这个图呢？我们通过一个简单实例来说明：
Train Data有890行。(请对应图中的上层部分）
每1次的fold，都会生成 712行小train， 178行小test。我们用Model 1来训练 712行的小train，然后预测 178行小test。预测的结果是长度为 178 的预测值。
这样的动作走5次！长度为178 的预测值 X 5 = 890 预测值，刚好和Train data长度吻合。这个890预测值是Model 1产生的，我们先存着，因为，一会让它将是第二层模型的训练来源。
重点：这一步产生的预测值我们可以转成 890 X 1 （890 行，1列），记作 P1 (大写P)
接着说 Test Data 有 418 行。(请对应图中的下层部分，对对对，绿绿的那些框框）
每1次的fold，712行小train训练出来的Model 1要去预测我们全部的Test Data（全部！因为Test Data没有加入5-fold，所以每次都是全部！）。此时，Model 1的预测结果是长度为418的预测值。
这样的动作走5次！我们可以得到一个 5 X 418 的预测值矩阵。然后我们根据行来就平均值，最后得到一个 1 X 418 的平均预测值。
重点：这一步产生的预测值我们可以转成 418 X 1 （418行，1列），记作 p1 (小写p)
走到这里，你的第一层的Model 1完成了它的使命。
第一层还会有其他Model的，比如Model 2，同样的走一遍，我们有可以得到 890 X 1 (P2) 和 418 X 1 (p2) 列预测值。
这样吧，假设你第一层有3个模型，这样你就会得到：
来自5-fold的预测值矩阵 890 X 3，（P1，P2， P3）和来自Test Data预测值矩阵 418 X 3，（p1, p2, p3）。
到第二层了
来自5-fold的预测值矩阵 890 X 3 作为你的Train Data，训练第二层的模型
来自Test Data预测值矩阵 418 X 3 就是你的Test Data，用训练好的模型来预测他们吧。

（1）堆叠平均模型类

class StackingAveragedModels(BaseEstimator, RegressorMixin, TransformerMixin):
    def __init__(self, base_models, meta_model, n_folds=5):
        self.base_models = base_models
        self.meta_model = meta_model
        self.n_folds = n_folds
   
    # We again fit the data on clones of the original models
    def fit(self, X, y):
        self.base_models_ = [list() for x in self.base_models]
        self.meta_model_ = clone(self.meta_model)
        kfold = KFold(n_splits=self.n_folds, shuffle=True, random_state=156)
        
        # Train cloned base models then create out-of-fold predictions
        # that are needed to train the cloned meta-model
        out_of_fold_predictions = np.zeros((X.shape[0], len(self.base_models)))
        for i, model in enumerate(self.base_models):
            for train_index, holdout_index in kfold.split(X, y):
                instance = clone(model)
                self.base_models_[i].append(instance)
                instance.fit(X[train_index], y[train_index])
                y_pred = instance.predict(X[holdout_index])
                out_of_fold_predictions[holdout_index, i] = y_pred
                
        # Now train the cloned  meta-model using the out-of-fold predictions as new feature
        self.meta_model_.fit(out_of_fold_predictions, y)
        return self
   
    #Do the predictions of all base models on the test data and use the averaged predictions as 
    #meta-features for the final prediction which is done by the meta-model
    def predict(self, X):
        meta_features = np.column_stack([
            np.column_stack([model.predict(X) for model in base_models]).mean(axis=1)
            for base_models in self.base_models_ ])
        return self.meta_model_.predict(meta_features)

class StackingAveragedModels(BaseEstimator, RegressorMixin, TransformerMixin):

def __init__(self, base_models, meta_model, n_folds=5):

self.base_models = base_models

self.meta_model = meta_model

self.n_folds = n_folds

# We again fit the data on clones of the original models

def fit(self, X, y):

self.base_models_ = [list() for x in self.base_models]

self.meta_model_ = clone(self.meta_model)

kfold = KFold(n_splits=self.n_folds, shuffle=True, random_state=156)

# Train cloned base models then create out-of-fold predictions

# that are needed to train the cloned meta-model

out_of_fold_predictions = np.zeros((X.shape[0], len(self.base_models)))

for i, model in enumerate(self.base_models):

for train_index, holdout_index in kfold.split(X, y):

instance = clone(model)

self.base_models_[i].append(instance)

instance.fit(X[train_index], y[train_index])

y_pred = instance.predict(X[holdout_index])

out_of_fold_predictions[holdout_index, i] = y_pred

# Now train the cloned meta-model using the out-of-fold predictions as new feature

self.meta_model_.fit(out_of_fold_predictions, y)

return self

#Do the predictions of all base models on the test data and use the averaged predictions as

#meta-features for the final prediction which is done by the meta-model

def predict(self, X):

meta_features = np.column_stack([

np.column_stack([model.predict(X) for model in base_models]).mean(axis=1)

for base_models in self.base_models_ ])

return self.meta_model_.predict(meta_features)

（2）堆叠平均模型得分
为了使两种方法具有可比性（通过使用相同数量的模型），我们只是平均Enet，KRR和Gboost，然后我们添加lasso作为元模型。

stacked_averaged_models = StackingAveragedModels(base_models = (ENet, GBoost, KRR),
                                                 meta_model = lasso)

score = rmsle_cv(stacked_averaged_models)
print("Stacking Averaged models score: {:.4f} ({:.4f})".format(score.mean(), score.std()))

stacked_averaged_models = StackingAveragedModels(base_models = (ENet, GBoost, KRR),

meta_model = lasso)

score = rmsle_cv(stacked_averaged_models)

print("Stacking Averaged models score: {:.4f} ({:.4f})".format(score.mean(), score.std()))

运行结果
Stacking Averaged models score: 0.0026 (0.0002)
由此可知，通过添加元学习器，我们再次获得更好的分数。
这个结果比简单求平均值的得分Averaged base models score: 0.0027 (0.0002)
更好！

（3）集成StackedRegressor，XGBoost和LightGBM
我们将XGBoost和LightGBM添加到之前定义的StackedRegressor中。
我们首先定义一个rmsle评估函数

def rmsle(y, y_pred):
    return np.sqrt(mean_squared_error(y, y_pred))

1 2	def rmsle(y, y_pred): return np.sqrt(mean_squared_error(y, y_pred))

（4）最终培训和预测

stacked_averaged_models.fit(train.values, y_train)
stacked_train_pred = stacked_averaged_models.predict(train.values)
stacked_pred = np.expm1(stacked_averaged_models.predict(test.values))
print(rmsle(y_train, stacked_train_pred))

stacked_averaged_models.fit(train.values, y_train)

stacked_train_pred = stacked_averaged_models.predict(train.values)

stacked_pred = np.expm1(stacked_averaged_models.predict(test.values))

print(rmsle(y_train, stacked_train_pred))

运行结果
0.00165651517208
（5）计算xgboost

model_xgb.fit(train, y_train)
xgb_train_pred = model_xgb.predict(train)
xgb_pred = np.expm1(model_xgb.predict(test))
print(rmsle(y_train, xgb_train_pred))

model_xgb.fit(train, y_train)

xgb_train_pred = model_xgb.predict(train)

xgb_pred = np.expm1(model_xgb.predict(test))

print(rmsle(y_train, xgb_train_pred))

运行结果
0.00860202560688
（6）计算LightGBM:

model_lgb.fit(train, y_train)
lgb_train_pred = model_lgb.predict(train)

lgb_pred = np.expm1(model_lgb.predict(test.values))
print(rmsle(y_train, lgb_train_pred))

model_lgb.fit(train, y_train)

lgb_train_pred = model_lgb.predict(train)

lgb_pred = np.expm1(model_lgb.predict(test.values))

print(rmsle(y_train, lgb_train_pred))

运行结果
0.0016023894136

'''RMSE on the entire Train data when averaging'''

print('RMSLE score on train data:')
print(rmsle(y_train,stacked_train_pred*0.70 +
               xgb_train_pred*0.15 + lgb_train_pred*0.15 ))

'''RMSE on the entire Train data when averaging'''

print('RMSLE score on train data:')

print(rmsle(y_train,stacked_train_pred*0.70 +

xgb_train_pred*0.15 + lgb_train_pred*0.15 ))

运行结果
RMSLE score on train data:
0.002334422609
由此看出，这似乎最简单的堆叠（或集成）方法确实提高了分数。
大家可参考：https://www.kaggle.com/serigne/stacked-regressions-top-4-on-leaderboard
集成学习非常成功，该算法不仅在挑战性的数据集上频频打破性能方面的记录，而且是 Kaggle 数据科学竞赛的获奖者常用的方法之一。
有关集成学习方法可参考：
https://zhuanlan.zhihu.com/p/25836678

17 降维简介

当特征选择完成后，可以直接训练模型了，但是可能由于特征矩阵过大，导致计算量大，训练时间长的问题，因此降低特征矩阵维度也是必不可少的。常见的降维方法除了以上提到的基于L1惩罚项的模型以外，另外还有主成分分析法（PCA）和线性判别分析（LDA），线性判别分析本身也是一个分类模型。PCA和LDA有很多的相似点，其本质是要将原始的样本映射到维度更低的样本空间中，但是PCA和LDA的映射目标不一样：PCA是为了让映射后的样本具有最大的发散性；而LDA是为了让映射后的样本有最好的分类性能。所以说PCA是一种无监督的降维方法，而LDA是一种有监督的降维方法。
PCA、LDA降维一般假设数据集为线性可分，如果用这两种方法，对线性不可分的数据集进行降维，效果往往不理想。本质上PCA和LDA还是一种线性变换。而线性不可分数据应该是很普遍的，对线性不可分数据集该如何进行降维呢？这里我们介绍一种核PCA方法，这样降维方法综合了核技巧及PCA思想，对非线性数据集降维有非常好的效果。
此外，这里我们还介绍SVD方法，这也是一种非常有效的降维方法。

17.1 PCA简介

主成分分析（Principal Components Analysis），简称PCA，是一种数据降维技术，用于数据预处理。一般我们获取的原始数据维度都很高，比如1000个特征，在这1000个特征中可能包含了很多无用的信息或者噪声，真正有用的特征才50个或更少，那么我们可以运用PCA算法将1000个特征降到50个特征。这样不仅可以去除无用的噪声，还能减少很大的计算量。
PCA算法是如何实现的？
简单来说，就是将数据从原特征空间转换到新的特征空间中，例如原始的空间是三维的(x,y,z)，x、y、z分别是原始空间的三个基，我们可以通过某种方法，用新的坐标系(a,b,c)来表示原始的数据，那么a、b、c就是新的基，它们组成新的特征空间。在新的特征空间中，可能所有的数据在c上的投影都接近于0，即可以忽略，那么我们就可以直接用(a,b)来表示数据，这样数据就从三维的(x,y,z)降到了二维的(a,b)。
问题是如何求新的基(a,b,c)?
一般步骤是这样的：
1）对原始数据集做标准化处理。
2）求协方差矩阵。
3）计算协方差矩阵的特征值和特征向量。
4）选择前k个最大的特征向量，k小于原数据集维度。
5）通过前k个特征向量组成了新的特征空间，设为W。
6)通过矩阵W,把原数据转换到新的k维特征子空间。

17.2 PCA算法实现

这里以葡萄酒数据为例，数据集特征如下：

数据来源于：https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data
1）对原数据集做标准化处理
导入需要的库及数据

import pandas as pd

df_wine = pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data', header=None)

df_wine.columns = ['Class label', 'Alcohol', 'Malic acid', 'Ash', 
'Alcalinity of ash', 'Magnesium', 'Total phenols', 
'Flavanoids', 'Nonflavanoid phenols', 'Proanthocyanins', 
'Color intensity', 'Hue', 'OD280/OD315 of diluted wines', 'Proline']

df_wine.head()

import pandas as pd

df_wine = pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data', header=None)

df_wine.columns = ['Class label', 'Alcohol', 'Malic acid', 'Ash',

'Alcalinity of ash', 'Magnesium', 'Total phenols',

'Flavanoids', 'Nonflavanoid phenols', 'Proanthocyanins',

'Color intensity', 'Hue', 'OD280/OD315 of diluted wines', 'Proline']

df_wine.head()

部分内容：

为便于后续处理，把数据集分为训练集和测试集，划分比例为7:3

from sklearn.cross_validation import train_test_split

X, y = df_wine.iloc[:, 1:].values, df_wine.iloc[:, 0].values
X_train, X_test, y_train, y_test = \
        train_test_split(X, y, test_size=0.3, random_state=0)

from sklearn.cross_validation import train_test_split

X, y = df_wine.iloc[:, 1:].values, df_wine.iloc[:, 0].values

X_train, X_test, y_train, y_test = \

train_test_split(X, y, test_size=0.3, random_state=0)

对原数据进行标准化处理

from sklearn.preprocessing import StandardScaler

sc = StandardScaler()
X_train_std = sc.fit_transform(X_train)
X_test_std = sc.fit_transform(X_test)

from sklearn.preprocessing import StandardScaler

sc = StandardScaler()

X_train_std = sc.fit_transform(X_train)

X_test_std = sc.fit_transform(X_test)

2) 求协方差矩阵
这里使用numpy.cov函数，求标准化后数据的协方差矩阵

import numpy as np
cov_mat = np.cov(X_train_std.T)

1 2	import numpy as np cov_mat = np.cov(X_train_std.T)

3）计算协方差矩阵的特征值和特征向量
使用np.linalg.eig函数，求协方差的特征值和特征向量

eigen_vals, eigen_vecs = np.linalg.eig(cov_mat)
print('\nEigenvalues \n%s' % eigen_vals)

1 2	eigen_vals, eigen_vecs = np.linalg.eig(cov_mat) print('\nEigenvalues \n%s' % eigen_vals)

得到13个特征向量：
Eigenvalues
[ 4.8923083 2.46635032 1.42809973 1.01233462 0.84906459 0.60181514 0.52251546 0.08414846 0.33051429 0.29595018 0.16831254 0.21432212 0.2399553 ]
要实现降维，我们可以选择前k个最多信息（或方差最大）特征向量组成新的子集，由于特征值的大小决定了特征向量的重要性，因此，可以通过对特征值的排序，获取前k个特征值。特征值λ_i的方差贡献率是指特征值λ_i与所有特征值和的比例：

我们可以通过numpy.cumsum函数计算累计方差。

tot = sum(eigen_vals)
var_exp = [(i / tot) for i in sorted(eigen_vals, reverse=True)]
cum_var_exp = np.cumsum(var_exp)

#然后用matplotlib各主成分的方差贡献率图形。

import matplotlib.pyplot as plt
%matplotlib inline
import matplotlib.font_manager as fm
myfont = fm.FontProperties(fname='/home/hadoop/anaconda3/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/simhei.ttf')

plt.bar(range(1, 14), var_exp, alpha=0.5, align='center',
        label='individual explained variance')
plt.step(range(1, 14), cum_var_exp, where='mid',
         label='cumulative explained variance')
plt.ylabel('方差贡献率',fontproperties=myfont,size=12)
plt.xlabel('主成分',fontproperties=myfont,size=12)
plt.legend(loc='best')
plt.tight_layout()
plt.show()

tot = sum(eigen_vals)

var_exp = [(i / tot) for i in sorted(eigen_vals, reverse=True)]

cum_var_exp = np.cumsum(var_exp)

#然后用matplotlib各主成分的方差贡献率图形。

import matplotlib.pyplot as plt

%matplotlib inline

import matplotlib.font_manager as fm

myfont = fm.FontProperties(fname='/home/hadoop/anaconda3/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/simhei.ttf')

plt.bar(range(1, 14), var_exp, alpha=0.5, align='center',

label='individual explained variance')

plt.step(range(1, 14), cum_var_exp, where='mid',

label='cumulative explained variance')

plt.ylabel('方差贡献率',fontproperties=myfont,size=12)

plt.xlabel('主成分',fontproperties=myfont,size=12)

plt.legend(loc='best')

plt.tight_layout()

plt.show()

从这个图可以看出第一个主成分占了方差总和的40%左右，前两个主成分占了近60%。
4）选择前k个最大的特征向量，k小于原数据集维度
首先，按特征值按降序排序

# 构成一个元组 (eigenvalue, eigenvector) 
eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:,i]) for i in range(len(eigen_vals))]

# Sort the (eigenvalue, eigenvector) tuples from high to low
eigen_pairs.sort(reverse=True)

# 构成一个元组 (eigenvalue, eigenvector)

eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:,i]) for i in range(len(eigen_vals))]

# Sort the (eigenvalue, eigenvector) tuples from high to low

eigen_pairs.sort(reverse=True)

5）通过前k个特征向量组成了新的特征空间，设为W。
为便于数据可视化，这里我们取k=2，实际上前2个特征值已占了总方差的近60%。

w = np.hstack((eigen_pairs[0][1][:, np.newaxis],
               eigen_pairs[1][1][:, np.newaxis]))
print('Matrix W:\n', w)

w = np.hstack((eigen_pairs[0][1][:, np.newaxis],

eigen_pairs[1][1][:, np.newaxis]))

print('Matrix W:\n', w)

这样我们就可得到一个由这两个特征向量构成的13*2矩阵W:
Matrix W:
[[ 0.14669811 0.50417079]
[-0.24224554 0.24216889]
[-0.02993442 0.28698484]
[-0.25519002 -0.06468718]
[ 0.12079772 0.22995385]
[ 0.38934455 0.09363991]
[ 0.42326486 0.01088622]
[-0.30634956 0.01870216]
[ 0.30572219 0.03040352]
[-0.09869191 0.54527081]
[ 0.30032535 -0.27924322]
[ 0.36821154 -0.174365 ]
[ 0.29259713 0.36315461]]

6）通过矩阵W,把原数据转换到新的k维特征子空间
通过这个特征矩阵W,把原样本x转换到PCA的子空间上，得到一个新样本x^，。
x^，=xW
训练集与W点积后，把这个训练集转换到包括两个主成分的子空间上。然后，把子空间的数据可视化。

X_train_pca = X_train_std.dot(w)
colors = ['r', 'b', 'g']
markers = ['s', 'x', 'o']

for l, c, m in zip(np.unique(y_train), colors, markers):
    plt.scatter(X_train_pca[y_train==l, 0], 
                X_train_pca[y_train==l, 1], 
                c=c, label=l, marker=m)

plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend(loc='lower left')
plt.tight_layout()
plt.show()

X_train_pca = X_train_std.dot(w)

colors = ['r', 'b', 'g']

markers = ['s', 'x', 'o']

for l, c, m in zip(np.unique(y_train), colors, markers):

plt.scatter(X_train_pca[y_train==l, 0],

X_train_pca[y_train==l, 1],

c=c, label=l, marker=m)

plt.xlabel('PC 1')

plt.ylabel('PC 2')

plt.legend(loc='lower left')

plt.tight_layout()

plt.show()

从以上图形可以看出，大部分数据沿PC1方向分布，而且可以线性划分，在可视化图形时，为便于标识点，这里采用了y_train标签信息。

我们用来6步来实现PCA，这个过程还是比较麻烦的，是否有更简单的方法呢？
有的，接下来我们介绍利用Scikit-learn中PCA类进行降维。

17.3 利用Scikit-learn进行主成分分析

我们将使用Scikit-learn中PCA对数据集进行预测处理，然后使用逻辑斯谛回归对转换后的数据进行分类，最后对数据进行可视化。
1）数据预处理

from sklearn.decomposition import PCA

pca = PCA()
X_train_pca = pca.fit_transform(X_train_std)
pca.explained_variance_ratio_

from sklearn.decomposition import PCA

pca = PCA()

X_train_pca = pca.fit_transform(X_train_std)

pca.explained_variance_ratio_

得到主成分数据：
array([ 0.37329648, 0.18818926, 0.10896791, 0.07724389, 0.06478595, 0.04592014, 0.03986936, 0.02521914, 0.02258181, 0.01830924, 0.01635336, 0.01284271, 0.00642076])

2）可视化主成分方差贡献率图

plt.bar(range(1, 14), pca.explained_variance_ratio_, alpha=0.5, align='center')
plt.step(range(1, 14), np.cumsum(pca.explained_variance_ratio_), where='mid')
plt.ylabel('Explained variance ratio')
plt.xlabel('Principal components')
plt.show()

plt.bar(range(1, 14), pca.explained_variance_ratio_, alpha=0.5, align='center')

plt.step(range(1, 14), np.cumsum(pca.explained_variance_ratio_), where='mid')

plt.ylabel('Explained variance ratio')

plt.xlabel('Principal components')

plt.show()

3）获取前2个主成分

pca = PCA(n_components=2)
X_train_pca = pca.fit_transform(X_train_std)
X_test_pca = pca.transform(X_test_std)
4）把训练集映射到主成分空间上，并可视化。
plt.scatter(X_train_pca[:,0], X_train_pca[:,1])
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.show()

pca = PCA(n_components=2)

X_train_pca = pca.fit_transform(X_train_std)

X_test_pca = pca.transform(X_test_std)

4）把训练集映射到主成分空间上，并可视化。

plt.scatter(X_train_pca[:,0], X_train_pca[:,1])

plt.xlabel('PC 1')

plt.ylabel('PC 2')

plt.show()

5）利用回归模型对数据进行分类。

from sklearn.linear_model import LogisticRegression

lr = LogisticRegression()
lr = lr.fit(X_train_pca, y_train)

from sklearn.linear_model import LogisticRegression

lr = LogisticRegression()

lr = lr.fit(X_train_pca, y_train)

6）为了更好看到分类后情况，这里我们定义一个函数plot_decision_regions，通过这个函数对决策区域数据可视化。

from matplotlib.colors import ListedColormap
def plot_decision_regions(X, y, classifier, resolution=0.02):

    # setup marker generator and color map
    markers = ('s', 'x', 'o', '^', 'v')
    colors = ('red', 'blue', 'lightgreen', 'gray', 'cyan')
    cmap = ListedColormap(colors[:len(np.unique(y))])

    # plot the decision surface
    x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, resolution),
                         np.arange(x2_min, x2_max, resolution))
    Z = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)
    Z = Z.reshape(xx1.shape)
    plt.contourf(xx1, xx2, Z, alpha=0.4, cmap=cmap)
    plt.xlim(xx1.min(), xx1.max())
    plt.ylim(xx2.min(), xx2.max())

    # plot class samples
    for idx, cl in enumerate(np.unique(y)):
        plt.scatter(x=X[y == cl, 0], y=X[y == cl, 1],
                    alpha=0.8, c=cmap(idx),
                    marker=markers[idx], label=cl)

from matplotlib.colors import ListedColormap

def plot_decision_regions(X, y, classifier, resolution=0.02):

# setup marker generator and color map

markers = ('s', 'x', 'o', '^', 'v')

colors = ('red', 'blue', 'lightgreen', 'gray', 'cyan')

cmap = ListedColormap(colors[:len(np.unique(y))])

# plot the decision surface

x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1

x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1

xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, resolution),

np.arange(x2_min, x2_max, resolution))

Z = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)

Z = Z.reshape(xx1.shape)

plt.contourf(xx1, xx2, Z, alpha=0.4, cmap=cmap)

plt.xlim(xx1.min(), xx1.max())

plt.ylim(xx2.min(), xx2.max())

# plot class samples

for idx, cl in enumerate(np.unique(y)):

plt.scatter(x=X[y == cl, 0], y=X[y == cl, 1],

alpha=0.8, c=cmap(idx),

marker=markers[idx], label=cl)

7)把训练数据转换到前两个主成分轴后生成决策区域图形

plot_decision_regions(X_train_pca, y_train, classifier=lr)
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend(loc='lower left')
plt.tight_layout()
# plt.savefig('./figures/pca3.png', dpi=300)
plt.show()

plot_decision_regions(X_train_pca, y_train, classifier=lr)

plt.xlabel('PC 1')

plt.ylabel('PC 2')

plt.legend(loc='lower left')

plt.tight_layout()

# plt.savefig('./figures/pca3.png', dpi=300)

plt.show()

对高维数据集进行降维除了PCA方法，还有线性判别分析（Linear Discriminant Analysis， LDA）、决策树、核主成分分析、SVD等等。

17.4 LDA 降维

LDA的基本概念与PCA类似，PCA是在数据集中找到方差最大的正交的主成分分量的轴。而LDA的目标是发现可以最优化分类的特征子空间。两者都是可以用于降维的线性转换方法，其中，PCA是无监督算法，LDA是监督算法。与PCA相比，LDA是一种更优越的用于分类的特征提取技术。
LDA的主要步骤：
（1）对d维数据集进行标准化处理（d为特征数量）
（2）对每一类别，计算d维的均值向量
（3）构造类间的散布矩阵S_B以及类内的散布矩阵S_W
（4）计算矩阵〖S_W〗^(-1) S_B的特征值所对应的特征向量，
（5）选取前k个特征值对应的特征向量，构造一个d x k维的转换矩阵W，其中特征向量以列的形式排列
（6）使用转换矩阵W将样本映射到新的特征子空间上.
以下还是以下葡萄酒数据为例，用代码实现以上各步：
（1）对d维数据集进行标准化处理

from sklearn.preprocessing import StandardScaler

X, y = df_wine.iloc[:, 1:].values, df_wine.iloc[:, 0].values
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0)
#对特征进行标准化处理
sc = StandardScaler()
X_train_std = sc.fit_transform(X_train)
X_test_std = sc.transform(X_test)

from sklearn.preprocessing import StandardScaler

X, y = df_wine.iloc[:, 1:].values, df_wine.iloc[:, 0].values

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0)

#对特征进行标准化处理

sc = StandardScaler()

X_train_std = sc.fit_transform(X_train)

X_test_std = sc.transform(X_test)

（2）对每一类别，计算d维的均值向量

#设置精度
np.set_printoptions(precision=4)

#求各类的平均值
mean_vecs = []
for label in range(1,4):
    mean_vecs.append(np.mean(X_train_std[y_train==label], axis=0))
    print('MV %s: %s\n' %(label, mean_vecs[label-1]))

#设置精度

np.set_printoptions(precision=4)

#求各类的平均值

mean_vecs = []

for label in range(1,4):

mean_vecs.append(np.mean(X_train_std[y_train==label], axis=0))

print('MV %s: %s\n' %(label, mean_vecs[label-1]))

运行结果
MV 1: [ 0.9259 -0.3091 0.2592 -0.7989 0.3039 0.9608 1.0515 -0.6306 0.5354 0.2209 0.4855 0.798 1.2017]

MV 2: [-0.8727 -0.3854 -0.4437 0.2481 -0.2409 -0.1059 0.0187 -0.0164 0.1095 -0.8796 0.4392 0.2776 -0.7016]

MV 3: [ 0.1637 0.8929 0.3249 0.5658 -0.01 -0.9499 -1.228 0.7436 -0.7652 0.979 -1.1698 -1.3007 -0.3912]
（3）构造类间的散布矩阵S_B以及类内的散布矩阵S_W
通过均值向量计算类内散布矩阵Sw：

通过累加各类别i的散布矩阵Si来计算：

d = 13 # number of features
S_W = np.zeros((d, d))
for label,mv in zip(range(1, 4), mean_vecs):
    class_scatter = np.zeros((d, d)) # scatter matrix for each class
    for row in X_train_std[y_train == label]:
        row, mv = row.reshape(d, 1), mv.reshape(d, 1) # make column vectors
        class_scatter += (row-mv).dot((row-mv).T)
    S_W += class_scatter                             # sum class scatter matrices

print('Within-class scatter matrix: %sx%s' % (S_W.shape[0], S_W.shape[1]))

d = 13 # number of features

S_W = np.zeros((d, d))

for label,mv in zip(range(1, 4), mean_vecs):

class_scatter = np.zeros((d, d)) # scatter matrix for each class

for row in X_train_std[y_train == label]:

row, mv = row.reshape(d, 1), mv.reshape(d, 1) # make column vectors

class_scatter += (row-mv).dot((row-mv).T)

S_W += class_scatter # sum class scatter matrices

print('Within-class scatter matrix: %sx%s' % (S_W.shape[0], S_W.shape[1]))

运行结果
Within-class scatter matrix: 13x13

计算各类标样本数

print('Class label distribution: %s' % np.bincount(y_train)[1:])

1	print('Class label distribution: %s' % np.bincount(y_train)[1:])

运行结果为：
Class label distribution: [40 49 35]
由此看出，各类记录数不很均匀，为此，需要对SB进行归一化处理：

d = 13 # number of features
S_W = np.zeros((d, d))
for label,mv in zip(range(1, 4), mean_vecs):
    class_scatter = np.cov(X_train_std[y_train==label].T)
    S_W += class_scatter
print('Scaled within-class scatter matrix: %sx%s' % (S_W.shape[0], S_W.shape[1]))

d = 13 # number of features

S_W = np.zeros((d, d))

for label,mv in zip(range(1, 4), mean_vecs):

class_scatter = np.cov(X_train_std[y_train==label].T)

S_W += class_scatter

print('Scaled within-class scatter matrix: %sx%s' % (S_W.shape[0], S_W.shape[1]))

运行结果
Scaled within-class scatter matrix: 13x13

计算类间散布矩阵：

mean_overall = np.mean(X_train_std, axis=0)
d = 13 # number of features
S_B = np.zeros((d, d))
for i,mean_vec in enumerate(mean_vecs):
    n = X_train[y_train==i+1, :].shape[0]
    mean_vec = mean_vec.reshape(d, 1) # make column vector
    mean_overall = mean_overall.reshape(d, 1) # make column vector
    S_B += n * (mean_vec - mean_overall).dot((mean_vec - mean_overall).T)

print('Between-class scatter matrix: %sx%s' % (S_B.shape[0], S_B.shape[1]))

mean_overall = np.mean(X_train_std, axis=0)

d = 13 # number of features

S_B = np.zeros((d, d))

for i,mean_vec in enumerate(mean_vecs):

n = X_train[y_train==i+1, :].shape[0]

mean_vec = mean_vec.reshape(d, 1) # make column vector

mean_overall = mean_overall.reshape(d, 1) # make column vector

S_B += n * (mean_vec - mean_overall).dot((mean_vec - mean_overall).T)

print('Between-class scatter matrix: %sx%s' % (S_B.shape[0], S_B.shape[1]))

运行结果
Between-class scatter matrix: 13x13

eigen_vals, eigen_vecs = np.linalg.eig(np.linalg.inv(S_W).dot(S_B))

1	eigen_vals, eigen_vecs = np.linalg.eig(np.linalg.inv(S_W).dot(S_B))

（5）选取前k个特征值对应的特征向量，构造一个d x k维的转换矩阵W，其中特征向量以列的形式排列

求得广义特征值之后，按照降序对特征值排序

# 生成特征值与特征向量构成的元组
eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:,i]) for i in range(len(eigen_vals))]

# Sort the (eigenvalue, eigenvector) tuples from high to low
eigen_pairs = sorted(eigen_pairs, key=lambda k: k[0], reverse=True)

# Visually confirm that the list is correctly sorted by decreasing eigenvalues

print('Eigenvalues in decreasing order:\n')
for eigen_val in eigen_pairs:
    print(eigen_val[0])

# 生成特征值与特征向量构成的元组

eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:,i]) for i in range(len(eigen_vals))]

# Sort the (eigenvalue, eigenvector) tuples from high to low

eigen_pairs = sorted(eigen_pairs, key=lambda k: k[0], reverse=True)

# Visually confirm that the list is correctly sorted by decreasing eigenvalues

print('Eigenvalues in decreasing order:\n')

for eigen_val in eigen_pairs:

print(eigen_val[0])

运行结果
Eigenvalues in decreasing order:

452.721581245
156.43636122
7.05575044266e-14
5.68434188608e-14
3.41129233161e-14
3.40797229523e-14
3.40797229523e-14
1.16775565372e-14
1.16775565372e-14
8.59477909861e-15
8.59477909861e-15
4.24523361436e-15
2.6858909629e-15
d x d维协方差矩阵的秩最大为d-1，得到两个非0的特征值。
与PCA一样，我们可视化各特征贡献率

tot = sum(eigen_vals.real)
discr = [(i / tot) for i in sorted(eigen_vals.real, reverse=True)]
cum_discr = np.cumsum(discr)

plt.bar(range(1, 14), discr, alpha=0.5, align='center',
        label='individual "discriminability"')
plt.step(range(1, 14), cum_discr, where='mid',
         label='cumulative "discriminability"')
plt.ylabel('"discriminability" ratio')
plt.xlabel('Linear Discriminants')
plt.ylim([-0.1, 1.1])
plt.legend(loc='best')
plt.tight_layout()
# plt.savefig('./figures/lda1.png', dpi=300)
plt.show()

tot = sum(eigen_vals.real)

discr = [(i / tot) for i in sorted(eigen_vals.real, reverse=True)]

cum_discr = np.cumsum(discr)

plt.bar(range(1, 14), discr, alpha=0.5, align='center',

label='individual "discriminability"')

plt.step(range(1, 14), cum_discr, where='mid',

label='cumulative "discriminability"')

plt.ylabel('"discriminability" ratio')

plt.xlabel('Linear Discriminants')

plt.ylim([-0.1, 1.1])

plt.legend(loc='best')

plt.tight_layout()

# plt.savefig('./figures/lda1.png', dpi=300)

plt.show()

运行结果

（6）使用转换矩阵W将样本映射到新的特征子空间上.
由上面两个新得到两个特征构成一个新矩阵

w = np.hstack((eigen_pairs[0][1][:, np.newaxis].real,
                      eigen_pairs[1][1][:, np.newaxis].real))
print('Matrix W:\n', w)

w = np.hstack((eigen_pairs[0][1][:, np.newaxis].real,

eigen_pairs[1][1][:, np.newaxis].real))

print('Matrix W:\n', w)

d x d维协方差矩阵的秩最大为d-1，得到两个非0的特征值。Matrix W:
[[-0.0662 -0.3797]
[ 0.0386 -0.2206]
[-0.0217 -0.3816]
[ 0.184 0.3018]
[-0.0034 0.0141]
[ 0.2326 0.0234]
[-0.7747 0.1869]
[-0.0811 0.0696]
[ 0.0875 0.1796]
[ 0.185 -0.284 ]
[-0.066 0.2349]
[-0.3805 0.073 ]
[-0.3285 -0.5971]]
将样本映射到新的特征空间

X_train_lda = X_train_std.dot(w)
colors = ['r', 'b', 'g']
markers = ['s', 'x', 'o']

for l, c, m in zip(np.unique(y_train), colors, markers):
    plt.scatter(X_train_lda[y_train==l, 0] * (-1), 
                X_train_lda[y_train==l, 1] * (-1), 
                c=c, label=l, marker=m)

plt.xlabel('LD 1')
plt.ylabel('LD 2')
plt.legend(loc='lower right')
plt.tight_layout()
# plt.savefig('./figures/lda2.png', dpi=300)
plt.show()

X_train_lda = X_train_std.dot(w)

colors = ['r', 'b', 'g']

markers = ['s', 'x', 'o']

for l, c, m in zip(np.unique(y_train), colors, markers):

plt.scatter(X_train_lda[y_train==l, 0] * (-1),

X_train_lda[y_train==l, 1] * (-1),

c=c, label=l, marker=m)

plt.xlabel('LD 1')

plt.ylabel('LD 2')

plt.legend(loc='lower right')

plt.tight_layout()

# plt.savefig('./figures/lda2.png', dpi=300)

plt.show()

运行结果

17.5 利用Scikit-learn进行LDA分析

下面我们利用scikit-learn中对LDA类的实现
这里先定义一个函数，后面需要用到

from matplotlib.colors import ListedColormap

def plot_decision_regions(X, y, classifier, resolution=0.02):

    # setup marker generator and color map
    markers = ('s', 'x', 'o', '^', 'v')
    colors = ('red', 'blue', 'lightgreen', 'gray', 'cyan')
    cmap = ListedColormap(colors[:len(np.unique(y))])

    # plot the decision surface
    x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, resolution),
                         np.arange(x2_min, x2_max, resolution))
    Z = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)
    Z = Z.reshape(xx1.shape)
    plt.contourf(xx1, xx2, Z, alpha=0.4, cmap=cmap)
    plt.xlim(xx1.min(), xx1.max())
    plt.ylim(xx2.min(), xx2.max())

    # plot class samples
    for idx, cl in enumerate(np.unique(y)):
        plt.scatter(x=X[y == cl, 0], y=X[y == cl, 1],
                    alpha=0.8, c=cmap(idx),
                    marker=markers[idx], label=cl)

from matplotlib.colors import ListedColormap

def plot_decision_regions(X, y, classifier, resolution=0.02):

# setup marker generator and color map

markers = ('s', 'x', 'o', '^', 'v')

colors = ('red', 'blue', 'lightgreen', 'gray', 'cyan')

cmap = ListedColormap(colors[:len(np.unique(y))])

# plot the decision surface

x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1

x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1

xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, resolution),

np.arange(x2_min, x2_max, resolution))

Z = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)

Z = Z.reshape(xx1.shape)

plt.contourf(xx1, xx2, Z, alpha=0.4, cmap=cmap)

plt.xlim(xx1.min(), xx1.max())

plt.ylim(xx2.min(), xx2.max())

# plot class samples

for idx, cl in enumerate(np.unique(y)):

plt.scatter(x=X[y == cl, 0], y=X[y == cl, 1],

alpha=0.8, c=cmap(idx),

marker=markers[idx], label=cl)

对数据先LDA处理，然后用逻辑回归进行分类。

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA

 
lda = LDA(n_components=2)
X_train_lda = lda.fit_transform(X_train_std, y_train)
# 逻辑回归在相对低维数据上的表现
from sklearn.linear_model import LogisticRegression
lr = LogisticRegression()
lr = lr.fit(X_train_lda, y_train)

plot_decision_regions(X_train_lda, y_train, classifier=lr)
plt.xlabel('LD 1')
plt.ylabel('LD 2')
plt.legend(loc='lower left')
plt.tight_layout()
# plt.savefig('./images/lda3.png', dpi=300)
plt.show()

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA

lda = LDA(n_components=2)

X_train_lda = lda.fit_transform(X_train_std, y_train)

# 逻辑回归在相对低维数据上的表现

from sklearn.linear_model import LogisticRegression

lr = LogisticRegression()

lr = lr.fit(X_train_lda, y_train)

plot_decision_regions(X_train_lda, y_train, classifier=lr)

plt.xlabel('LD 1')

plt.ylabel('LD 2')

plt.legend(loc='lower left')

plt.tight_layout()

# plt.savefig('./images/lda3.png', dpi=300)

plt.show()

运行结果

还有几个点划分错误，下面通过正则化，效果将更好

X_test_lda = lda.transform(X_test_std)

plot_decision_regions(X_test_lda, y_test, classifier=lr)
plt.xlabel('LD 1')
plt.ylabel('LD 2')
plt.legend(loc='lower left')
plt.tight_layout()
# plt.savefig('./images/lda4.png', dpi=300)
plt.show()

X_test_lda = lda.transform(X_test_std)

plot_decision_regions(X_test_lda, y_test, classifier=lr)

plt.xlabel('LD 1')

plt.ylabel('LD 2')

plt.legend(loc='lower left')

plt.tight_layout()

# plt.savefig('./images/lda4.png', dpi=300)

plt.show()

运行结果

17.6使用核PCA降维

前面我们介绍了两种降维方法，PCA及LDA.这两种方法，如果用于线性不可分数据集上进行分类，效果往往不很理想，原因是通过他们无法把线性不可分数据集变为线性可分数据集。如果遇到线性不可分数据集（这样的数据集往往比较普遍），有什么好方法，既降维，又可把线性不可分数据集变为线性可分数据集？
在SVM中，我们了解到核函数的神奇，把可以通过把线性不可分的数据集映射到一个高维空间，变得线性可分。基于这点，如果我们在降维时也采用核技术是否也可以呢？可以的，这就是接下来我们要介绍的内容---核PCA.
核PCA=核技术+PCA,具体步骤如下：
（1）计算核矩阵，也就是计算任意两个训练样本。这里以向基核函数（RBF）为例
经向基函数核（又称高斯核）为：

得到以下矩阵：

（2）对核矩阵K进行中心化处理

其中,是n*n的矩阵，n=训练集样本数，中每个元素都等于.l_n中的每个元素都是1/n
（3）求核矩阵的特征向量，并按降序排列，提取前k个特征向量。
不同于标准PCA，这里的特征向量并不是主成分轴。
下面我们根据以上三个步骤，实现一个核PCA。借助SciPy和NumPy，其实实现核PCA很简单。

from scipy.spatial.distance import pdist, squareform
from scipy import exp
from scipy.linalg import eigh
import numpy as np

def rbf_kernel_pca(X, gamma, n_components):
    """
    RBF kernel PCA implementation.

    Parameters
    ------------
    X: {NumPy ndarray}, shape = [n_samples, n_features]
        
    gamma: float
      Tuning parameter of the RBF kernel
        
    n_components: int
      Number of principal components to return

    Returns
    ------------
     X_pc: {NumPy ndarray}, shape = [n_samples, k_features]
       Projected dataset   

    """
    # Calculate pairwise squared Euclidean distances
    # in the MxN dimensional dataset.
    sq_dists = pdist(X, 'sqeuclidean')

    # Convert pairwise distances into a square matrix.
    mat_sq_dists = squareform(sq_dists)

    # Compute the symmetric kernel matrix.
    K = exp(-gamma * mat_sq_dists)

    # Center the kernel matrix.
    N = K.shape[0]
    one_n = np.ones((N,N)) / N
    K = K - one_n.dot(K) - K.dot(one_n) + one_n.dot(K).dot(one_n)

    # Obtaining eigenpairs from the centered kernel matrix
    # numpy.eigh returns them in sorted order
    eigvals, eigvecs = eigh(K)

    # Collect the top k eigenvectors (projected samples)
    X_pc = np.column_stack((eigvecs[:, -i]
                            for i in range(1, n_components + 1)))

    return X_pc

from scipy.spatial.distance import pdist, squareform

from scipy import exp

from scipy.linalg import eigh

import numpy as np

def rbf_kernel_pca(X, gamma, n_components):

"""

RBF kernel PCA implementation.

Parameters

------------

X: {NumPy ndarray}, shape = [n_samples, n_features]

gamma: float

Tuning parameter of the RBF kernel

n_components: int

Number of principal components to return

Returns

------------

X_pc: {NumPy ndarray}, shape = [n_samples, k_features]

Projected dataset

"""

# Calculate pairwise squared Euclidean distances

# in the MxN dimensional dataset.

sq_dists = pdist(X, 'sqeuclidean')

# Convert pairwise distances into a square matrix.

mat_sq_dists = squareform(sq_dists)

# Compute the symmetric kernel matrix.

K = exp(-gamma * mat_sq_dists)

# Center the kernel matrix.

N = K.shape[0]

one_n = np.ones((N,N)) / N

K = K - one_n.dot(K) - K.dot(one_n) + one_n.dot(K).dot(one_n)

# Obtaining eigenpairs from the centered kernel matrix

# numpy.eigh returns them in sorted order

eigvals, eigvecs = eigh(K)

# Collect the top k eigenvectors (projected samples)

X_pc = np.column_stack((eigvecs[:, -i]

for i in range(1, n_components + 1)))

return X_pc

下面以一分离同心数据集为例，分别用PCA和核PCA对数据集进行处理，然后处理后的结果，具体请看以下代码及生成的图形：

from sklearn.datasets import make_circles

X, y = make_circles(n_samples=1000, random_state=123, noise=0.1, factor=0.2)

plt.scatter(X[y==0, 0], X[y==0, 1], color='red', marker='^', alpha=0.5)
plt.scatter(X[y==1, 0], X[y==1, 1], color='blue', marker='o', alpha=0.5)

plt.tight_layout()
plt.show()

from sklearn.datasets import make_circles

X, y = make_circles(n_samples=1000, random_state=123, noise=0.1, factor=0.2)

plt.scatter(X[y==0, 0], X[y==0, 1], color='red', marker='^', alpha=0.5)

plt.scatter(X[y==1, 0], X[y==1, 1], color='blue', marker='o', alpha=0.5)

plt.tight_layout()

plt.show()

这是典型线性不可数据集，现在我们分别用PCA及核PCA进行处理。
（1）用PCA处理，然后进行分类

from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt

scikit_pca = PCA(n_components=2)
X_spca = scikit_pca.fit_transform(X)

plt.figure( figsize=(5,3))

plt.scatter(X_spca[y==0, 0], X_spca[y==0, 1], 
            color='red', marker='^', alpha=0.5)
plt.scatter(X_spca[y==1, 0], X_spca[y==1, 1],
            color='blue', marker='o', alpha=0.5)


plt.xlabel('PC1')
plt.ylabel('PC2')

plt.show()

from sklearn.decomposition import PCA

from sklearn.preprocessing import StandardScaler

import matplotlib.pyplot as plt

scikit_pca = PCA(n_components=2)

X_spca = scikit_pca.fit_transform(X)

plt.figure( figsize=(5,3))

plt.scatter(X_spca[y==0, 0], X_spca[y==0, 1],

color='red', marker='^', alpha=0.5)

plt.scatter(X_spca[y==1, 0], X_spca[y==1, 1],

color='blue', marker='o', alpha=0.5)

plt.xlabel('PC1')

plt.ylabel('PC2')

plt.show()

（2）用核PCA处理，然后进行分类

X_kpca = rbf_kernel_pca(X, gamma=15, n_components=2)

plt.figure( figsize=(5,3))
plt.scatter(X_kpca[y==0, 0], X_kpca[y==0, 1], 
            color='red', marker='^', alpha=0.5)
plt.scatter(X_kpca[y==1, 0], X_kpca[y==1, 1],
            color='blue', marker='o', alpha=0.5)

plt.xlabel('PC1')
plt.ylabel('PC2')

plt.show()

X_kpca = rbf_kernel_pca(X, gamma=15, n_components=2)

plt.figure( figsize=(5,3))

plt.scatter(X_kpca[y==0, 0], X_kpca[y==0, 1],

color='red', marker='^', alpha=0.5)

plt.scatter(X_kpca[y==1, 0], X_kpca[y==1, 1],

color='blue', marker='o', alpha=0.5)

plt.xlabel('PC1')

plt.ylabel('PC2')

plt.show()

（3）使用sklearn实现核PCA
源数据的图形为

这里通过核PCA把该数据变为线性可分数据集，实现代码如下：

from sklearn.decomposition import KernelPCA
from sklearn.datasets import make_moons

X, y = make_moons(n_samples=100, random_state=123)
scikit_kpca = KernelPCA(n_components=2, kernel='rbf', gamma=15)
X_skernpca = scikit_kpca.fit_transform(X)

plt.scatter(X_skernpca[y==0, 0], X_skernpca[y==0, 1], 
            color='red', marker='^', alpha=0.5)
plt.scatter(X_skernpca[y==1, 0], X_skernpca[y==1, 1], 
            color='blue', marker='o', alpha=0.5)

plt.xlabel('PC1')
plt.ylabel('PC2')
plt.tight_layout()
plt.show()

from sklearn.decomposition import KernelPCA

from sklearn.datasets import make_moons

X, y = make_moons(n_samples=100, random_state=123)

scikit_kpca = KernelPCA(n_components=2, kernel='rbf', gamma=15)

X_skernpca = scikit_kpca.fit_transform(X)

plt.scatter(X_skernpca[y==0, 0], X_skernpca[y==0, 1],

color='red', marker='^', alpha=0.5)

plt.scatter(X_skernpca[y==1, 0], X_skernpca[y==1, 1],

color='blue', marker='o', alpha=0.5)

plt.xlabel('PC1')

plt.ylabel('PC2')

plt.tight_layout()

plt.show()

17.7 SVD矩阵分解

（1）SVD奇异值分解的定义
假设有一个mxn矩阵，如果存在一个分解

其中U为的mxm酉矩阵，∑为mxn的半正定对角矩阵，除了对角元素不为0，其他元素都为0，并且对角元素是从大到小排列的，前面的元素比较大，后面的很多元素接近0。这些对角元素就是奇异值。V^T为V的共轭转置矩阵，且为nxn的酉矩阵。这样的分解称为的奇异值分解，对角线上的元素称为奇异值，U称为左奇异矩阵，V^T称为右奇异矩阵。
SVD在信息检索（隐性语义索引）、图像压缩、推荐系统、金融等领域都有应用。
（2）SVD奇异值分解与特征值分解的关系
特征值分解与SVD奇异值分解的目的都是提取一个矩阵最重要的特征。然而，特征值分解只适用于方阵，而SVD奇异值分解适用于任意的矩阵，不一定是方阵。

这里M^T M和MM^T都是方阵，UU^T和VV^T都是单位矩阵，V是M^T M的特征向量，U是MM^T的特征向量。
（3）SVD奇异值分解的作用和意义
奇异值分解最大的作用就是数据的降维，当然，还有其他很多的作用，这里主要讨论数据的降维，对于mxn的M矩阵，进行奇异值分解

取其前k个非零奇异值，可以还原原来的矩阵，即前k个非零奇异值对应的奇异向量代表了矩阵的主要特征。可以表示为

17.8 用Python实现SVD，并用于图像压缩
（1）首先读取一张图片（128*128*3）：

#!python
#  -*- coding:utf-8 -*-
from PIL import Image
import os
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
 
if __name__ == '__main__':
    mpl.rcParams['font.sans-serif'] = [u'simHei']
    mpl.rcParams['axes.unicode_minus'] = False
    A = Image.open('02.jpg')
    a = np.array(A)  #转换成矩阵

#!python

# -*- coding:utf-8 -*-

from PIL import Image

import os

import numpy as np

import matplotlib as mpl

import matplotlib.pyplot as plt

if __name__ == '__main__':

mpl.rcParams['font.sans-serif'] = [u'simHei']

mpl.rcParams['axes.unicode_minus'] = False

A = Image.open('02.jpg')

a = np.array(A) #转换成矩阵

（2）然后可以利用python的numpy库对彩色图像的3个通道进行SVD分解

numpy库中有SVD分解函数：np.linalg.svd
#由于是彩色图像，所以3通道。a的最内层数组为三个数，分别表示RGB，用来表示一个像素
u_r, sigma_r, v_r = np.linalg.svd(a[:, :, 0])
u_g, sigma_g, v_g = np.linalg.svd(a[:, :, 1])
u_b, sigma_b, v_b = np.linalg.svd(a[:, :, 2])

numpy库中有SVD分解函数：np.linalg.svd

#由于是彩色图像，所以3通道。a的最内层数组为三个数，分别表示RGB，用来表示一个像素

u_r, sigma_r, v_r = np.linalg.svd(a[:, :, 0])

u_g, sigma_g, v_g = np.linalg.svd(a[:, :, 1])

u_b, sigma_b, v_b = np.linalg.svd(a[:, :, 2])

（3）然后便可以根据需要压缩图像（丢弃分解出来的三个矩阵中的数据），利用的奇异值个数越少，则压缩的越厉害。下面来看一下不同程度压缩后，重构图像的清晰度：

plt.figure(facecolor = 'w', figsize = (10, 10))
# 奇异值个数依次取：1,2,...,16。来看看一下效果
K = 16
for k in range(1, K + 1):
    R = restore(u_r, sigma_r, v_r, k)
    G = restore(u_g, sigma_g, v_g, k)
    B = restore(u_b, sigma_b, v_b, k)
    I = np.stack((R, G, B), axis = 2)
    # 将图片重构后的显示出来
    plt.subplot(4, 4, k)
    plt.imshow(I)
    plt.axis('off')
    plt.title(u'奇异值个数：%d' %  k)

plt.suptitle(u'SVD与图像分解', fontsize = 20)
plt.tight_layout(0.1, rect = (0, 0, 1, 0.92))
plt.show()

plt.figure(facecolor = 'w', figsize = (10, 10))

# 奇异值个数依次取：1,2,...,16。来看看一下效果

K = 16

for k in range(1, K + 1):

R = restore(u_r, sigma_r, v_r, k)

G = restore(u_g, sigma_g, v_g, k)

B = restore(u_b, sigma_b, v_b, k)

I = np.stack((R, G, B), axis = 2)

# 将图片重构后的显示出来

plt.subplot(4, 4, k)

plt.imshow(I)

plt.axis('off')

plt.title(u'奇异值个数：%d' % k)

plt.suptitle(u'SVD与图像分解', fontsize = 20)

plt.tight_layout(0.1, rect = (0, 0, 1, 0.92))

plt.show()

（4）其中restore函数定义为

def restore(u, sigma, v, k):
    m = len(u)
    n = len(v)
    a = np.zeros((m, n))
    # 重构图像
    a = np.dot(u[:, :k], np.diag(sigma[:k])).dot(v[:k, :])
    # 上述语句等价于：
    # for i in range(k):
    #     ui = u[:, i].reshape(m, 1)
    #     vi = v[i].reshape(1, n)
    #     a += sigma[i] * np.dot(ui, vi)
    a[a < 0] = 0 a[a > 255] = 255
    return np.rint(a).astype('uint8')

def restore(u, sigma, v, k):

m = len(u)

n = len(v)

a = np.zeros((m, n))

# 重构图像

a = np.dot(u[:, :k], np.diag(sigma[:k])).dot(v[:k, :])

# 上述语句等价于：

# for i in range(k):

# ui = u[:, i].reshape(m, 1)

# vi = v[i].reshape(1, n)

# a += sigma[i] * np.dot(ui, vi)

a[a < 0] = 0 a[a > 255] = 255

return np.rint(a).astype('uint8')

Python技术交流与分享

分享技术平台

月度归档：2018年09月

18.1集成学习概述

18.1.1 Bagging

18.1.2 boosting

18.1.3 Stacking

18.1.3 .1Stacking中元分类层

18.1.3.2Stacking的几种方法

18.2投票分类器（VotingClassifier）

18.2.1多数投票分类（MajorityVote Class）

18.2.1.1Iris数据集概述

18.2.2多数投票分类（MajorityVote Class）

18.3自适应分类器（Adaboost）

18.4 Xgboost简介

18.4.1简介

18.4.2 xgboost 实例

18.4.2.1数据探索及数据预处理

18.4.2.2创建模型

18.4.2.3看各模型的性能

18.4.2.4堆叠（或集成）各模型

18.4.2.5添加元模型

17 降维简介

17.1 PCA简介

17.2 PCA算法实现

17.3 利用Scikit-learn进行主成分分析

17.4 LDA 降维

17.5 利用Scikit-learn进行LDA分析

17.6使用核PCA降维

17.7 SVD矩阵分解