Python 数据挖掘

AI 培训练习的笔记

下载数据

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

import os
import urllib.request

# https://archive.ics.uci.edu/ml/datasets/Iris
if not os.path.exists('iris.csv'):
    urlString = 'http://aima.cs.berkeley.edu/data/iris.csv'
    irisFile = urllib.request.urlopen(urlString)
    localIrisFile = open('iris.csv', 'wb')
    localIrisFile.write(irisFile.read())
    localIrisFile.close()

我们使用了 urllib 类库获取伯克利大学网站的一个数据文件，并使把它保存到本地磁盘。

数据包含鸢尾花（iris）数据集，这是一个包含了三种鸢尾花（山鸢尾 setosa、维吉尼亚鸢尾 virginica 和变色鸢尾 versicolor）的各 50 个数据样本的多元数据集。

每个样本都有四个特征（或者说变量），即花萼（sepal）和花瓣（petal）的长度和宽度，单位为厘米。

数据集以 CSV（逗号分割值）的格式存储。CSV 文件可以很方便的转化并把其中的信息存储为适合的数据结构。此数据集有 5 列。格式如下：

5.0,3.3,1.4,0.2,setosa
7.0,3.2,4.7,1.4,versicolor
7.1,3.0,5.9,2.1,virginica

数据导入

from numpy import genfromtxt

# 读取前四列特征 花萼（sepal）的长度和宽度花瓣（petal）的长度和宽度
data = genfromtxt('iris.csv', delimiter=',', usecols=(0, 1, 2, 3))
# 读取最后一列标注
target = genfromtxt('iris.csv', delimiter=',', usecols=(4,), dtype=str)

print(data.shape)
# (150, 4)
print(target.shape)
# (150,)
print(set(target))
# {'versicolor', 'virginica', 'setosa'}

上面的代码中创建了一个包含特征值的矩阵和一个包含样本类型的向量。我们可以通过查看我们加载的数据结构的 shape 来确认数据集的大小，也可以查看我们有多少种样本类型以及它们的名字。

当我们处理新数据的时候，一项很重要的任务是尝试去理解数据包含的信息以及它的组织结构。可视化可以灵活生动的展示数据，帮助我们深入理解数据。

可视化

from pylab import plot, show

plot(data[target == 'setosa', 0], data[target == 'setosa', 2], 'bo')
plot(data[target == 'versicolor', 0], data[target == 'versicolor', 2], 'r+')
plot(data[target == 'virginica', 0], data[target == 'virginica', 2], 'g*')
show()

from pylab import figure, subplot, hist, xlim

xmin = min(data[:, 0])
xmax = max(data[:, 0])
figure()  # 新图层
subplot(411)  # 图层分区
hist(data[target == 'setosa', 0], color='b', alpha=.7)  # 绘制直方图
xlim(xmin, xmax)  # 设置坐标系
subplot(412)
hist(data[target == 'versicolor', 0], color='r', alpha=.7)
xlim(xmin, xmax)
subplot(413)
hist(data[target == 'virginica', 0], color='g', alpha=.7)
xlim(xmin, xmax)
subplot(414)
hist(data[:, 0], color='y', alpha=.7)
xlim(xmin, xmax)
show()

图 1 中有 150 个点，不同的颜色代表不同的类型；蓝色点代表山鸢尾，红色点代表变色鸢尾，绿色点代表维吉尼亚鸢尾。

另一种常用的查看数据的方法是分特性绘制直方图。在图 2 中，数据被分为三类，我们可以比较每一类的分布特征。第二段代码绘制数据中每一类型的第一个特性（花萼的长度）

分类

from numpy import zeros

# 将分类用数字代替
t = zeros(len(target))
t[target == 'setosa'] = 1
t[target == 'versicolor'] = 2
t[target == 'virginica'] = 3

# 使用连续型朴素贝叶斯分类
from sklearn.naive_bayes import GaussianNB

classifier = GaussianNB()
classifier.fit(data, t)

print("预测：", (classifier.predict([data[0]])))
# 预测： [ 1.]
print("结果：", t[0])
# 结果： 1.0

分类是一个数据挖掘方法，用于把一个数据集中的样本数据分配给各个目标类。实现这个方法的模块叫做分类器。使用分类器需要以下两步：训练和分类。训练是指采集已知其特定类归属的数据并基于这些数据创建分类器。 分类是指使用通过这些已知数据建立的分类器来处理未知的数据，以判断未知数据的分类情况。

Sklearn 类库包含很多分类器的实现，第一段代码使用高斯朴素贝叶斯来分析我们在最开始载入的鸢尾花数据，包含山鸢尾、变色鸢尾和维吉尼亚鸢尾。最后把字符串数组转型成整型数据

分类器可以由 predict 方法完成，并且只要输出一个样例就可以很简单的检测。第二段代码中 predicted 类包含了一个正确的样本（山鸢尾），但是在广泛的样本上评估分类器并且使用非训练环节的数据测试是很重要的。最终我们通过从源数据集中随机抽取样本把数据分为训练集和测试集。我们将会使用训练集的数据来训练分类器，并使用测试集的数据来测试分类器。train_test_split 方法正是实现此功能的。

from sklearn import model_selection

train, test, t_train, t_test = model_selection.train_test_split(data, t, test_size=0.4, random_state=0)

classifier.fit(train, t_train)  # 训练
print(classifier.score(test, t_test))  # 测试
# 0.933333333333

数据集被分一分为二，测试集被指定为源数据的40%（命名为 test_size），我们用它反复训练我们的分类器并输出精确度。

在此例中，我们的精确度为 93%。一个分类器的精确度是通过正确分类样本的数量除以总样本的数量得出的。也就是说，它意味着我们正确预测的比例。

from sklearn.metrics import confusion_matrix

print(confusion_matrix(classifier.predict(test), t_test))
# [[16  0  0]
#  [ 0 23  4]
#  [ 0  0 17]]

另一个估计分类器表现的工具叫做混淆矩阵。在此矩阵中每列代表一个预测类的实例，每行代表一个实际类的实例。使用它可以很容易的计算和打印矩阵

在这个混淆矩阵中我们可以看到所有山鸢尾和维吉尼亚鸢尾都被正确的分类了，但是实际上应该是 26 个的变色鸢尾，系统却预测其中三个是维吉尼亚鸢尾。如果我们牢记所有正确的猜测都在表格的对角线上，那么观测表格的错误就很容易了，即对角线以外的非零值。

from sklearn.metrics import classification_report

print(classification_report(classifier.predict(test), t_test, target_names=['setosa', 'versicolor', 'virginica']))
#              precision    recall  f1-score   support
#
#      setosa       1.00      1.00      1.00        16
#  versicolor       1.00      0.85      0.92        27
#   virginica       0.81      1.00      0.89        17
#
# avg / total       0.95      0.93      0.93        60

Precision：正确预测的比例，属于该类的被准确分类为该类的比例。100% 也有可能包含别的样本被错误的归类到本类目下。

Recall（或者叫真阳性率）：正确识别的比例，被划分到此类下面的是否确实属于此类。100% 也有可能漏掉属于该类却没有划分到该类的样本。

F1-Score：precision和recall的调和平均数

以上仅仅只是给出用于支撑测试分类的数据量。当然，分割数据、减少用于训练的样本数以及评估结果等操作都依赖于配对的训练集和测试集的随机选择。如果要切实评估一个分类器并与其它的分类器作比较的话，我们需要使用一个更加精确的评估模型，例如 Cross Validation。该模型背后的思想很简单：多次将数据分为不同的训练集和测试集，最终分类器评估选取多次预测的平均值。这次，sklearn 为我们提供了运行模型的方法

from sklearn.model_selection import cross_val_score

scores = cross_val_score(classifier, data, t, cv=6)
print(scores)
# [ 0.92592593  1.          0.91666667  0.91666667  0.95833333  1.        ]

from numpy import mean

print(mean(scores))
# 0.952932098765

如上所见，输出是每次模型迭代产生的精确度的数组。我们可以很容易计算出平均精确度

聚类

通常我们的数据上不会有标签告诉我们它的样本类型；我们需要分析数据，把数据按照它们的相似度标准分成不同的群组，群组（或者群集）指的是相似样本的集合。这种分析被称为无监督数据分析。最著名的聚类工具之一叫做 k-means 算法，如下所示：

划分:K-means,k-MEDOIDs

层次聚类：CURE，BIRCH

密度聚类：DBSCAN

from sklearn.cluster import KMeans

kmeans = KMeans(n_clusters=3, init='random')
c = kmeans.fit_predict(data)

from sklearn.metrics import completeness_score, homogeneity_score

print('完整性得分', completeness_score(t, c))
# 完整性得分 0.764986151449
print('同质性得分', homogeneity_score(t, c))
# 同质性得分 0.751485402199

上述片段运行 k-measn 算法并把数据分为三个群集（参数 k 所指定的）。现在我们可以使用模型把每一个样本分配到三个群集中。我们可以估计群集的结果，与使用完整性得分和同质性得分计算而得的标签作比较

同质性 homogeneity：每个群集只包含单个类的成员；

完整性 completeness：给定类的所有成员都分配给同一个群集。

figure()
subplot(211)
plot(data[t == 1, 0], data[t == 1, 2], 'bo')
plot(data[t == 2, 0], data[t == 2, 2], 'ro')
plot(data[t == 3, 0], data[t == 3, 2], 'go')
subplot(212)
plot(data[c == 0, 0], data[c == 0, 2], 'bo', alpha=.7)
plot(data[c == 1, 0], data[c == 1, 2], 'ro', alpha=.7)
plot(data[c == 2, 0], data[c == 2, 2], 'go', alpha=.7)
show()

我们可以集群可视化并和带有真实标签的做可视化比较。

观察此图我们可以看到，底部左侧的群集可以被 k-means 完全识别，然而顶部的两个群集有部分识别错误。

回归

回归是一个用于预测变量之间函数关系的方法。例如，我们有两个变量，一个被认为是解释，一个被认为是依赖。我们希望使用模型描述两者的关系。当这种关系是一条线的时候就称为线性回归。

为了应用线性回归我们建立一个由上所述的综合数据集：

from numpy.random import rand

x = rand(40, 1)  # 解释变量
y = x * x * x + rand(40, 1) / 5  # 依赖变量

使用在 sklear.linear_model 模块中的 LinearRegression 模型。该模型可以通过计算每个数据点到拟合线的垂直差的平方和，找到平方和最小的最佳拟合线。使用方法和我们之前遇到的实现 sklearn 的模型类似。然后通过把拟合线和实际数据点画在同一幅图上来评估结果：

from sklearn.linear_model import LinearRegression

linreg = LinearRegression()
linreg.fit(x, y)

from numpy import linspace, matrix

figure()
xx = linspace(0, 1, 40)
plot(x, y, 'o', xx, linreg.predict(matrix(xx).T), '--r')
show()

观察该图我们可以得出结论：拟合线从数据点中心穿过，并可以确定是增长的趋势。

使用均方误差来量化模型和原始数据的拟合度：

from sklearn.metrics import mean_squared_error

print('拟合度：', mean_squared_error(linreg.predict(x), y))
# 拟合度： 0.0175373207578

相关

我们通过研究相关性来理解成对的变量之间是否相关，相关性的强弱。此类分析帮助我们精确定位被依赖的重要变量。最好的相关方法是皮尔逊积矩相关系数。它是由两个变量的协方差除以它们的标准差的乘积计算而来。我们将鸢尾花数据集的变量两两组合计算出其系数如下所示：

from numpy import corrcoef

corr = corrcoef(data.T)  # .T gives the transpose
print(corr)
# [[ 1.         -0.10936925  0.87175416  0.81795363]
#  [-0.10936925  1.         -0.4205161  -0.35654409]
#  [ 0.87175416 -0.4205161   1.          0.9627571 ]
#  [ 0.81795363 -0.35654409  0.9627571   1.        ]]

corrcoef 方法通过输入行为变量列为观察值的矩阵，计算返回相关系数的对称矩阵。该矩阵的每个元素代表着两个变量的相关性。当值一起增长时相关性为正。当一个值减少而另一个只增加时相关性为负。特别说明，1 代表完美的正相关，0 代表不相关，-1 代表完美的负相关。

当变量数增长时我们可以使用伪彩色点很方便的可视化相关矩阵：

from pylab import pcolor, colorbar, xticks, yticks
from numpy import arange

figure()
pcolor(corr)
colorbar()  # add
# 在坐标系上标识名称
xticks(arange(0.5, 4.5), ['sepal length', 'sepal width', 'petal length', 'petal width'], rotation=-20)
yticks(arange(0.5, 4.5), ['sepal length', 'sepal width', 'petal length', 'petal width'], rotation=-20)
show()

看图右侧的彩条，我们可以把颜色点关联到数值上。在本例中，红色被关联为最高的正相关，我们可以看出我们数据集的最强相关是“花瓣宽度”和“花瓣长度”这两个变量。

降维

单独使用该方法，我们只能看到数据集的部分数据视图。既然我们可以同时绘制的最高维度数为3，将整个数据集嵌入一系列维度并建立一个整体可视化视图是很有必要的。这个嵌入过程就被称作降维。最著名的降维技术之一就是主成分分析（PCA）。该技术把数据变量转换为等量或更少的不相关变量，称为主成分（PCs）。

们实例化了一个PCA对象，用于计算前两个主成分。然后如往常一样绘制结果：

from sklearn.decomposition import PCA

figure()
pca = PCA(n_components=2)
pcad = pca.fit_transform(data)
plot(pcad[target == 'setosa', 0], pcad[target == 'setosa', 1], 'bo')
plot(pcad[target == 'versicolor', 0], pcad[target == 'versicolor', 1], 'ro')
plot(pcad[target == 'virginica', 0], pcad[target == 'virginica', 1], 'go')
show()

可以注意到上图和第一章提到的有些相似，不过这次变色鸢尾（红色的）和维吉尼亚鸢尾（绿色的）的间隔更清晰了。

print(pca.explained_variance_ratio_)
print(1 - sum(pca.explained_variance_ratio_))
data_inv = pca.inverse_transform(pcad)
print(abs(sum(sum(data - data_inv))))
for i in range(1, 5):
    pca = PCA(n_components=i)
    pca.fit(data)
    print(sum(pca.explained_variance_ratio_) * 100, '%')

PCA将空间数据方差最大化，我们可以通过方差比判断 PCs 包含的信息量：pca.explained_variance_ratio_

现在我们知道第一个PC占原始数据的92%的信息量而第二个占剩下的5%。我们还可以输出在转化过程中丢失的信息量：pca.explained_variance_ratio_

在本例中我们损失了2%的信息量。此时，我们可以是应用逆变换还原原始数据：pca.inverse_transform

可以证明的是，由于信息丢失逆变换不能给出准确的原始数据。我们可以估算逆变换的结果和原始数据的相似度：abs(sum(sum(data - data_inv)))

可以看出原始数据和逆变换计算出的近似值之间的差异接近于零。通过改变主成分的数值来计算我们能够覆盖多少信息量是很有趣的。PCs 用得越多，信息覆盖就越全，不过这段分析有助于我们理解保存一段特定的信息需要哪些组件。例如，从上述片段可以看出，只要使用三个 PCs 就可以覆盖鸢尾花数据集的几乎 100% 的信息。

完整代码

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

# 下载数据
import os
import urllib.request

# https://archive.ics.uci.edu/ml/datasets/Iris
if not os.path.exists('iris.csv'):
    urlString = 'http://aima.cs.berkeley.edu/data/iris.csv'
    irisFile = urllib.request.urlopen(urlString)
    localIrisFile = open('iris.csv', 'wb')
    localIrisFile.write(irisFile.read())
    localIrisFile.close()

# 数据导入
from numpy import genfromtxt

# 读取前四列特征 花萼（sepal）的长度和宽度花瓣（petal）的长度和宽度
data = genfromtxt('iris.csv', delimiter=',', usecols=(0, 1, 2, 3))
# 读取最后一列标注
target = genfromtxt('iris.csv', delimiter=',', usecols=(4,), dtype=str)

print(data.shape)
# (150, 4)
print(target.shape)
# (150,)
print(set(target))
# {'versicolor', 'virginica', 'setosa'}

from pylab import plot, show

plot(data[target == 'setosa', 0], data[target == 'setosa', 2], 'bo')
plot(data[target == 'versicolor', 0], data[target == 'versicolor', 2], 'r+')
plot(data[target == 'virginica', 0], data[target == 'virginica', 2], 'g*')
# show()

from pylab import figure, subplot, hist, xlim

xmin = min(data[:, 0])
xmax = max(data[:, 0])
figure()  # 新图层
subplot(411)  # 图层分区
hist(data[target == 'setosa', 0], color='b', alpha=.7)  # 绘制直方图
xlim(xmin, xmax)  # 设置坐标系
subplot(412)
hist(data[target == 'versicolor', 0], color='r', alpha=.7)
xlim(xmin, xmax)
subplot(413)
hist(data[target == 'virginica', 0], color='g', alpha=.7)
xlim(xmin, xmax)
subplot(414)
hist(data[:, 0], color='y', alpha=.7)
xlim(xmin, xmax)
# show()

# 分类
from numpy import zeros

# 将分类用数字代替
t = zeros(len(target))
t[target == 'setosa'] = 1
t[target == 'versicolor'] = 2
t[target == 'virginica'] = 3

# 使用连续型朴素贝叶斯分类
from sklearn.naive_bayes import GaussianNB

classifier = GaussianNB()
classifier.fit(data, t)

print('预测：', (classifier.predict([data[0]])))
# 预测： [ 1.]
print('结果：', t[0])
# 结果： 1.0

from sklearn import model_selection

train, test, t_train, t_test = model_selection.train_test_split(data, t, test_size=0.4, random_state=0)

classifier.fit(train, t_train)  # 训练
print(classifier.score(test, t_test))  # 测试
# 0.933333333333

from sklearn.metrics import confusion_matrix

print(confusion_matrix(classifier.predict(test), t_test))
# [[16  0  0]
#  [ 0 23  4]
#  [ 0  0 17]]

from sklearn.metrics import classification_report

print(classification_report(classifier.predict(test), t_test, target_names=['setosa', 'versicolor', 'virginica']))
#              precision    recall  f1-score   support
#
#      setosa       1.00      1.00      1.00        16
#  versicolor       1.00      0.85      0.92        27
#   virginica       0.81      1.00      0.89        17
#
# avg / total       0.95      0.93      0.93        60

from sklearn.model_selection import cross_val_score

scores = cross_val_score(classifier, data, t, cv=6)
print(scores)
# [ 0.92592593  1.          0.91666667  0.91666667  0.95833333  1.        ]

from numpy import mean

print('平均精确度', mean(scores))
# 平均精确度 0.952932098765

# 聚类
from sklearn.cluster import KMeans

kmeans = KMeans(n_clusters=3, init='random')
c = kmeans.fit_predict(data)

from sklearn.metrics import completeness_score, homogeneity_score

print('完整性得分', completeness_score(t, c))
# 完整性得分 0.764986151449
print('同质性得分', homogeneity_score(t, c))
# 同质性得分 0.751485402199

figure()
subplot(211)
plot(data[t == 1, 0], data[t == 1, 2], 'bo')
plot(data[t == 2, 0], data[t == 2, 2], 'ro')
plot(data[t == 3, 0], data[t == 3, 2], 'go')
subplot(212)
plot(data[c == 0, 0], data[c == 0, 2], 'bo', alpha=.7)
plot(data[c == 1, 0], data[c == 1, 2], 'ro', alpha=.7)
plot(data[c == 2, 0], data[c == 2, 2], 'go', alpha=.7)
# show()

# 回归
from numpy.random import rand

x = rand(40, 1)  # 解释变量
y = x * x * x + rand(40, 1) / 5  # 依赖变量

from sklearn.linear_model import LinearRegression

linreg = LinearRegression()
linreg.fit(x, y)

from numpy import linspace, matrix

figure()
xx = linspace(0, 1, 40)
plot(x, y, 'o', xx, linreg.predict(matrix(xx).T), '--r')
show()

from sklearn.metrics import mean_squared_error

print('拟合度：', mean_squared_error(linreg.predict(x), y))
# 拟合度： 0.0175373207578

# 相关
from numpy import corrcoef

corr = corrcoef(data.T)  # .T gives the transpose
print(corr)
# [[ 1.         -0.10936925  0.87175416  0.81795363]
#  [-0.10936925  1.         -0.4205161  -0.35654409]
#  [ 0.87175416 -0.4205161   1.          0.9627571 ]
#  [ 0.81795363 -0.35654409  0.9627571   1.        ]]

from pylab import pcolor, colorbar, xticks, yticks
from numpy import arange

figure()
pcolor(corr)
colorbar()  # add
# 在坐标系上标识名称
xticks(arange(0.5, 4.5), ['sepal length', 'sepal width', 'petal length', 'petal width'], rotation=-20)
yticks(arange(0.5, 4.5), ['sepal length', 'sepal width', 'petal length', 'petal width'], rotation=-20)
# show()

# 降维
from sklearn.decomposition import PCA

figure()
pca = PCA(n_components=2)
pcad = pca.fit_transform(data)
plot(pcad[target == 'setosa', 0], pcad[target == 'setosa', 1], 'bo')
plot(pcad[target == 'versicolor', 0], pcad[target == 'versicolor', 1], 'ro')
plot(pcad[target == 'virginica', 0], pcad[target == 'virginica', 1], 'go')
show()

print(pca.explained_variance_ratio_)
print(1 - sum(pca.explained_variance_ratio_))
data_inv = pca.inverse_transform(pcad)
print(abs(sum(sum(data - data_inv))))
for i in range(1, 5):
    pca = PCA(n_components=i)
    pca.fit(data)
    print(sum(pca.explained_variance_ratio_) * 100, '%')