分类
程序和算法

搜索引擎中的Page Rank算法

PageRank算法

PageRank算法总的来说就是预先给每个网页一个PR值(下面用PR值指代PageRank值),由于PR值物理意义上为一个网页被访问概率,所以一般是\(\frac{1}{N}\),其中N为网页总数。另外,一般情况下,所有网页的PR值的总和为1。如果不为1的话也不是不行,最后算出来的不同网页之间PR值的大小关系仍然是正确的,只是不能直接地反映概率了。

所以PageRank算法实际上就是预先给定PR值后,通过每个网页之间的链接关系不断迭代,直至达到平稳分布为止。

各个网页的PR值之间的关系一般情况下表示为如下的式子:

\(PR(p_i)=α\sum_{p_j∈M_{p_i}}\frac{PR(p_j)}{L(p_j)}+\frac{(1-α)}{N}\)

其中\(M_{p_i}\)是所有对\(p_i\)网页有出链的网页集合;\(L(p_j)\)是网页的出链数目;\(N\)是网页总数;\(α\)是阻尼系数,即用户离开当前网页重新输入网址访问的概率,一般取0.85。

根据这一关系不断迭代,当算法收敛的时候,得到的PR值即使每个网页的PR排序值。


算法实现

digraph类

由于一个搜索引擎中的网页数不胜数,使用一个矩阵存储各个网页之间的链接关系是不可能的,这里采用有向图类来优化稀疏矩阵的存储。

  • 类的方法一览表
名字 参数 说明
init_nodes self, nodes 初始化向图图的节点
init_edges self, edges 初始化向图图的边
add_edge self, edge 增加向图图中的边
neighbors self, node 返回参数节点的所有后继节点
incidents self, node 返回参数节点的所有前驱节点
  • digraph类源代码
class digraph:
 def init_nodes(self, nodes):
  self.nodes = nodes

 def init_edges(self, edges):
  self.edges = edges

 def add_edge(self, edge):
  self.edges.append(edge)

 def neighbors(self, node):
  neighbors = []
  for edge in self.edges:
   if edge[0] == node:
    neighbors.append(edge[1])
  return neighbors

 def incidents(self, node):
  incidents = []
  for edge in self.edges:
   if edge[1] == node:
    incidents.append(edge[0])
  return incidents

PRIterator类

PRIterator类是用来迭代计算PR值的类,由一个有向图和一个阻尼系数构造,阻尼系数即公式中的α。默认的最大迭代次数是100,确定迭代是否结束,即PR值是否收敛的参数为0.00001,即ϵ。

只有一个方法——page_rank方法,来迭代计算PR值直至收敛。同时如果一个节点没有任何的后继节点,那么相当于所有的节点都是它的后继节点,即一个网页没有出链的情况下,下一秒访问任何网页的概率都是相等的。

  • PRIterator类源代码
class PRIterator:
 __doc__ = '''计算一张图中的PR值'''

 def __init__(self, dg, damping_factor):
  self.damping_factor = damping_factor  # 阻尼系数,即α
  self.max_iterations = 100             # 最大迭代次数
  self.min_delta = 0.00001              # 确定迭代是否结束的参数,即ϵ
  self.graph = dg

 def page_rank(self):
  #  先将图中没有出链的节点改为对所有节点都有出链
  for node in self.graph.nodes:
   if len(self.graph.neighbors(node)) == 0:
    for node2 in self.graph.nodes:
     self.graph.add_edge((node, node2))

  nodes = self.graph.nodes
  graph_size = len(nodes)

  if graph_size == 0:
   return {}
  page_rank = dict.fromkeys(nodes, 1.0 / graph_size)  # 给每个节点赋予初始的PR值
  damping_value = (1.0 - self.damping_factor) / graph_size  # 公式中的(1−α)/N部分

  flag = False
  for i in range(self.max_iterations):
   change = 0
   for node in nodes:
    rank = 0
    for incident_page in self.graph.incidents(node):  # 遍历所有"入射"的页面
     rank += self.damping_factor * (page_rank[incident_page] / len(self.graph.neighbors(incident_page)))
    rank += damping_value
    change += abs(page_rank[node] - rank)  # 绝对值
    page_rank[node] = rank

   print("This is NO.%s iteration" % (i + 1))

   if change < self.min_delta:
    flag = True
    break
  if flag:
   print("finished in %s iterations!" % node)
  else:
   print("finished out of 100 iterations!")
  return page_rank

主函数

主函数的作用就是读取网页的链接数据,初始化有向图,迭代计算PR值直至收敛,最后输出计算结果。

测试的数据为维基百科网站上所有网页的链接数据,由于数据共有103689条,这里不再给出详细数据,文末附有数据集及源码,可执行文件下载链接。这里给出一部分样例的截图如下。数据格式为:网页ID——链接到的网页ID。

在主程序中,需要输入测试数据的路径和阻尼系数。运行结果如下图。

  • 主程序源代码
if __name__ == '__main__':
 # python NewPageRankDemo.py

 # 输入数据和输出结果所在的文件夹位置
 directory = input('输入数据和输出结果所在的文件夹位置:\n示例输入:C:\\Users\\37618\\Desktop\\PageRank\\\n')
 damping_factor = input('输入阻尼系数α(一般取0.85):\n')
 damping_factor = float(damping_factor)

 # 储存训练样本数据
 f = open(directory + 'WikiData.txt')

 nodes = []
 edges = []

 for v in f:
  tmp = v.split("\t")
  tmp[1] = tmp[1].split("\n")[0]
  if nodes.count(tmp[0]) == 0:
   nodes.append(tmp[0])
  if nodes.count(tmp[1]) == 0:
   nodes.append(tmp[1])
  edges.append((tmp[0], tmp[1]))

 # 输出预测结果
 import sys

 output = sys.stdout
 outputfile = open(directory + 'result(' +str(damping_factor)+ ').txt', 'w')
 sys.stdout = outputfile

 dg = digraph()

 dg.init_nodes(nodes)

 dg.init_edges(edges)

 pr = PRIterator(dg, damping_factor)
 page_ranks = pr.page_rank()

 print("The final page rank is\n")
 for page_rank in page_ranks:
  print(page_rank, ' ', page_ranks.get(page_rank))

 outputfile.close()
 sys.stdout = output

示例输出如下图。

可以看到,程序在31次迭代后在ID为8274的网页处收敛,并输出了所有网页的PR值。

[mdx_warning title=”特别注意”]10万条数据的迭代非常占内存,请配置垃圾的笔记本们慎重尝试![/mdx_warning]


源代码链接

源代码链接,附带测试数据集及测试输出,源代码及说明,可执行文件

分类
程序和算法

大数据和机器学习 基础篇 分类 支持向量机SVM

分类算法是机器学习中的一个重点,也是人们常说的“有监督的学习”。这是一种利用一系列已知类别的样本来对模型进行训练调整分类器的参数,使其达到所要求性能的过程,也成为监督训练或有教师学习。

注:本文中用到的Python及其模块安装教程参见


支持向量机SVM

支持向量机SVM是一种比较抽象的算法概念,全称是Support Vector Machine,它可以用来做模式识别,分类或者回归的机器学习。

前面介绍过机器学习是为了解决样本的具体分类映射的问题,构造一个算法,把已知样本的特征和分类情况做一个逻辑映射关系,这样碰到样本时就能用这个算法把它进行分类了。

例如,大于零的实数叫正数,小于零的实数叫负数:

\(
属性=
\left\{
\begin{array}{c}
Positive(x>0)\\
\\
Negative(x<0) \end{array} \right.
\)

但这个过程其实并不是一个机器学习过程,都是人来告诉计算机判定的定义如何,然后计算机根据这个判定的定义来处理每一个待定的对象。这里面计算机确实没有学习的过程。

下面我们把上面的例子改变一下:


年龄和好坏

假设某公司有任务,给客户分类,看看什么样的用户质量比较高;假设客户信息只有年龄这一项,客户信息表如下:

  • 表1 客户信息表

[mdx_table header=”true”]
客户编号
客户年龄
客户质量
客户编号
客户年龄
客户质量
—–
XXXX
34

XXXX
30

—–
XXXX
33

XXXX
25
不好
—–
XXXX
32

XXXX
23
不好
—–
XXXX
31

XXXX
22
不好
—–
XXXX
30

XXXX
18
不好
[/mdx_table]

  • 客户信息的数轴表示:

从上图可以看出,年龄在30以上的客户质量都好,在25一下的都不好。那么就可以考虑在30和25中间切一刀,一边是好一边是不好。也就是从27.5分开。

  • 以27.5为分界:

例如,来了一个27岁的客户,那么他应该属于客户质量好的客户。但是如果没有办法“一刀切”怎么办?例如客户信息如下表:

  • 表2 客户信息表

[mdx_table header=”true”]
客户编号
客户年龄
客户质量
客户编号
客户年龄
客户质量
—–
XXXX
34

XXXX
30

—–
XXXX
33

XXXX
25
不好
—–
XXXX
32

XXXX
23

—–
XXXX
31
不好
XXXX
22
不好
—–
XXXX
30

XXXX
18
不好
[/mdx_table]

  • 客户信息的数轴表示:

这个比较麻烦,因为没法做到“一刀切”。那么现在有以下两种选择:

  • 把这两类全部标出来。
\(
属性=
\left\{
\begin{array}{c}
Good(x≥30,x≠31)\\
\\
Good(x=23)\\
\\
NotGood(x≤25,x≠23)\\
\\
NotGood(x=31)
\end{array}
\right.
\)

这是一个分段函数,虽然准确但是啰嗦。在实际情况中一般可能会有10000个甚至更多的客户,所以这个函数写起来会相当麻烦。所以这种标记方法很容易产生过拟的问题。

  • 一刀切。反正只要一刀切就会非常简洁,如果一刀切下去虽然两边的类都不纯或者一边的类不纯,但是只要不纯的程度在能容忍的范围内就可以。
\(
属性=
\left\{
\begin{array}{c}
Good(x≥27.5)\\
\\
NotGood(x<27.5)\\ \end{array} \right.
\)

这次还是从27.5切,大于等于27.5就算“好”,小于27.5就算“不好”。那么在这个分类中“好”类的不纯度为\(\frac{1}{6}\);“不好”类的不纯度为\(\frac{1}{4}\)。如果说在不做这种数据分析的情况下,发展10个客户,有6个是“好”客户,4个是“不好”的客户;现在改进后虽然有一定的误判率,约为16.7%——由于分类不纯的问题,但是直接过滤掉了一定的对象,发展10个客户,有8.3个是“好”客户,1.7个是“不好”的客户。从数值上看,方案策略的提升还是有改进的。

总结一下,有以下几个关键点:

  • 关键点1:切下去的点在SVM算法体系里叫“超平面”。这个“超平面”是一个抽象的面概念,在一维空间里就是一个点,用\(x+A=0\)的点来表示;二维空间里就是一条线,用\(Ax+By+C=0\)的直线来表示;三维空间里就是一个面,用\(Ax+By+Cz+D=0\)的平面来表示;四维空间了就用\(Ax+By+Cz+Dα+E=0\)来表示;以此类推……上述4个方程都可以变形为:
    \( x=-A\)
    \( y=-\frac{A}{B}x-\frac{C}{B}\)
    \( z=-\frac{A}{C}x-\frac{B}{C}y-\frac{D}{C}\)
    \( α=-\frac{A}{D}x-\frac{B}{D}y-\frac{C}{D}z-\frac{E}{D}\)
  • 关键点2:过拟问题。一般来说,设计分类器都是要尽量避免过拟的。过拟会给归纳过程带来很大的麻烦,而且在应用的过程中也非常不方便,只要精确度达到标准就足够了。
  • 关键点3:不纯度问题。不纯度和精确度是一对矛盾,精确度越高那么不纯度就越低,反之,不纯度越高精确度就越低。分类器的研究和调整的过程就是一个精度和成本平衡的过程,所以并不是不纯度越低越好,而是在实际生产中操作成本一样的情况下,不纯度越低越好。

定义距离

在一个平面直角坐标系中,有一些样本点作为训练点,一些被标记为类别\(X\),一些被标记为非类别\(X\)。如何想把两种类别的点区分开,最朴素的想法是,如果能找到一条直线\(Ax+By+C=0\)能够恰好把它们区分成两个部分就最好了。

如何求解恰当的\(A,B,C\),恰当的标准是让类别\(X\)中与该直线最近的样本点的距离和非类别\(X\)中的样本点与该直线的距离最大。

在平面直角坐标系中,如果有一条直线\(Ax+By+C=0\),那么点\((x_0,y_0)\)到该直线的距离如下:

\( d=\frac{|Ax_0+By_0+C|}{\sqrt{A^2+B^2}}\)

与数轴上类似:“\(x-27.5>0\)的都是分类为‘好’的样本,\(x-27.5<0\)的则都不是分类为‘好’的样本”。如果以\(Ax+By+C=0\)这样一个方程做分隔,可以发现,所有\(X\)类别中的样本点都满足\(Ax+By+C>0\),而所有非\(X\)类别中的样本点都满足\(Ax+By+C<0\),这就是最开始构造这条直线的目的。 推广到N维空间,这个超平面的公式可以简写为: \( g(v)=wv+b\) 这里w是系数,v是样本向量,b是常数。 而N维空间上的距离为: \( d=\frac{|g(v_0)|}{||w||}\)


升维

如果样本是线性不可分的,那么就需要升维来解决这个问题——这才是SVM算法最为吸引人的部分。

先来看一个一维数据的例子,来说明升维的概念。

假设在数轴上给出一些数据,其中[-2,2]区间内的被标记为分类1,其余的都是分类0,我们可以把分段函数写成如下的形式:

\(
f(x)=
\left\{
\begin{array}{c}
1(-x^2+4>0)\\
\\
0(-x^2+4≤0)
\end{array}
\right.
\)

或者可以认为:

\(
f(y)=
\left\{
\begin{array}{c}
1(y>0)\\
\\
0(y≤0)
\end{array}
\right.\)

\(\left\{
\begin{array}{c}
y=-x^2+4\\
\\
z=f(y)
\end{array}
\right.
\)

这实际上就是\(y=-x^2+4\)这个函数在\(y=0(x轴)\)这条直线上的投影把样本点分开。

  • \(y=-x^2+4\)的图形:

在二维空间中也有类似的方式。

例如,样本向量v距离原点(0,0)的距离为1以内分类被标记为0,其余都为1。同样是线性不可分的,但是可以构造一个这样的函数:

\(
f(x,y)=
\left\{
\begin{array}{c}
1(x^2+y^2≥1)\\
\\
0(x^2+y^2<1) \end{array} \right.
\)

或者可以这么认为:

\(
f(z)=
\left\{
\begin{array}{c}
1(z≥1)\\
\\
0(z<1) \end{array} \right.
\)

\(
\left\{
\begin{array}{c}
z=x^2+y^2\\
\\
α=f(z)
\end{array}
\right.
\)

可以看到,在一维空间上解决线性不可分问题是把函数映射到二维空间,使得一维空间上的分类边界是二维空间上的分类函数在以为空间上的投影;而在二维空间上解决线性不可分问题是把函数映射到三维空间,使得二维空间上的分类边界是三维空间上的分类函数在二维空间上的投影。那么所有的n维空间上的线性不可分问题都可以考虑映射到n+1维上去构造分类函数,使得它在n维空间上的投影能够将两个类别分开。

这个构造的过程在SVM里是使用核函数(Kernel)来进行的。核函数的目的很单纯,即只要在当前维度的样本是线性不可分的,就一律映射到更高的维度上去,在更高的维度上找到超平面,得到超平面方程。函数本身表示的只是一个变量代换关系。


示例

下面还是用客户信息的列表给出示例。

  • 表3 客户信息表

[mdx_table header=”true”]
客户编号
客户年龄
客户质量
客户编号
客户年龄
客户质量
—–
XXXX
34

XXXX
30

—–
XXXX
33
不好
XXXX
25
不好
—–
XXXX
32

XXXX
23

—–
XXXX
31
不好
XXXX
22
不好
—–
XXXX
30

XXXX
18
不好
[/mdx_table]

在这个例子中,客户年龄和客户质量之间明显是没办法做线性分隔了。这是可以用SVM来做分类。
在Python的Scikit-learn库中,用到的是一个叫做SVC的类,SVC是Support Vector Classification的缩写,即支持分类向量机。SVC支持的核函数包括linear(线性核函数),poly(多项式核函数),rbf(径向核函数),sigmoid(神经元激活核函数),precomputed(自定义核函数),默认使用rbf(径向核函数)。

示例代码如下:

#支持向量机SVM
from sklearn import svm

#年龄
X=[[34],[33],[32],[31],[30],[30],[25],[23],[22],[18]]
#质量
y=[1,0,1,0,1,1,0,1,0,1]

#把训练数据和对应的分类放入分类器中进行训练
#这里使用rbf(linear,poly,rbf,sigmoid,precomputed)
clf=svm.SVC(kernel='rbf').fit(X,y)

#预测年龄29的客户的质量
p=[[29]]
print(clf.predict(p))

预测年龄为29的客户的质量,输出为1,表示客户质量为“好”。

在这个例子中,使用的是rbf核函数,这也是最适合做非线性关系分类标准的首选核函数。如果这几种核函数实在不知道该用哪个,那就在实际场景中多做对比测试,看看哪一种的正确率最高即可。


#小结

SVM解决问题的方法有以下几步:

  1. 把所有的样本和其对应的分类标记交给算法进行训练。
  2. 如果发现线性可分,那就直接找出超平面。
  3. 如果发现线性不可分,那就映射到n+1维空间,找出超平面。
  4. 最后得到超平面的表达式,也就是分类函数。

想了解更多关于大数据和机器学习:

分类
程序和算法

大数据和机器学习 基础篇 分类 隐马尔可夫模型

分类算法是机器学习中的一个重点,也是人们常说的“有监督的学习”。这是一种利用一系列已知类别的样本来对模型进行训练调整分类器的参数,使其达到所要求性能的过程,也成为监督训练或有教师学习。

注:本文中用到的Python及其模块安装教程参见


隐马尔可夫模型

隐马尔可夫模型(Hidden Markov Model,HMM)最初由L. E. Baum发表在20世纪70年代一系列的统计学论文中,随后在语言识别,自然语言处理以及生物信息等领域体现了很大的价值。

首先提出马尔可夫链的概念:

  • 在观察一个系统变化的时候,它下一个状态(第n+1个状态)如何的概率只需要观察和统计当前状态(第n个状态)即可正确得出。

如果在一个完整的观察过程中有一些状态的转换,如下图中的\(y_1\)到\(y_n\)。在观察中\(y_1\)到\(y_n\)的状态存在一个客观的转化规律,但是没办法直接观测到,观测到的是每个\(y\)状态下的输出\(x\),即\(x_1\)到\(x_n\)。需要通过\(x_1\)到\(x_n\)这些输出值来进行模型建立和计算状态转移概率。

  • 隐马尔可夫链示意图:

为了比较容易理解整个过程,下面举一个例子,假设有3中不同的骰子:
第一个骰子是常见的骰子(称这个骰子为D6),6个面,每个面(1,2,3,4,5,6)出现的概率是\(\frac{1}{6}\)。
第二个骰子是一个四面体(称这个骰子为D4),每个面(1,2,3,4)出现的概率是\(\frac{1}{4}\)。
第三个骰子有8个面(称这个骰子为D8),每个面(1,2,3,4,5,6,7,8)出现的概率是\(\frac{1}{8}\)。

  • 三种骰子和掷骰子可能产生的结果:

先随机选择一个骰子,然后再用它掷出一个数字,并记录下这个选择和数字。先从三个骰子里挑一个,挑到每个骰子的概率都是\(\frac{1}{3}\)。然后掷骰子,得到一个数字,1,2,3,4,5,6,7,8中的一个。不停地重复上述过程,会得到一串数字,每个数字都是1,2,3,4,5,6,7,8中的一个。

例如,投掷骰子10次,可能得到这么一串数字:1,6,3,5,2,7,3,5,2,4,这串数字叫做可见状态链,也就是记录的这组数字,可是前面介绍的\(x_n\)。但是在隐马尔可夫模型中,不仅有可见状态链,还有一串隐含状态链。在这个例子里,这串隐含状态链就是选出的骰子的序列。例如本次中,隐含状态链是:D6,D8,D8,D6,D4,D8,D6,D6,D4,D8,如下图所示:

  • 隐马尔可夫模型示意图:

一般来说,HMM中的马尔科夫链就是指隐含状态链,因为实际是隐含状态之间存在转换概率(Transition Probability)

其实对于HMM来说,如果提前知道所有隐含状态之间的转换概率和所有隐含状态到所有可见状态之间的输出概率,进行模拟是相当容易的。但应用HMM模型时,往往缺失一部分信息,有时候知道骰子有几种,每种骰子是什么,但是不知道掷出来的骰子序列;有时候只是看到了很多次掷骰子的结果,剩下的什么都不知道。如何应用算法去估计这些缺失的信息,就成了一个很重要的问题。

HMM模型相关的算法主要分为3类,分别解决3中问题:

  • 问题1:知道骰子有几种(隐含状态数量),每种骰子是什么(转换概率),根据掷骰子掷出的结果(可见状态链),想知道每次掷出来的都是哪种骰子(隐含状态链)。

这个问题在语音识别领域叫做解码问题。这个问题其实有两种解法,会给出两个不同的答案。每个答案都正确,只是这些答案的意义不一样。第一种解法求最大似然状态路径,通俗地说,就是求一串骰子序列,这串骰子序列产生观测结果的概率最大。第二种解法不是求一组骰子序列,而是求每次掷出的骰子分别是某种骰子的概率。

例如,看到结果后,可以求得第一次掷骰子是D4的概率是0.5,D6的概率是0.3,D8的概率是0.2。

  • 问题2:知道骰子有几种(隐含状态数量),每种骰子是什么(转换概率),根据掷骰子掷出的结果(可见状态链),想知道掷出这个结果的概率。

这个问题看似意义不大,因为掷出来的结果很多时候都对应了一个比较大的概率。这个问题的目的其实是检测观察到的结果和已知的模型是否吻合。如果很多次结果都对应了比较小的概率,那么说明已知的模型很有可能是错的,有人偷偷把骰子换了。

  • 问题3:知道骰子有几种(隐含状态数量),不知道每种骰子是什么(转换概率),观测到很多次掷骰子的结果(可见状态链),想反推出每种骰子是什么(转换概率)。这个问题很重要,因为这是最常见的情况。很多时候只有可见结果,不知道HMM模型中的参数,需要从可见结果估计出这些参数,这是建模的一个必要步骤。

维特比算法

维特比算法用于解决前面提到的问题1的解最大似然路径的问题。

还是考虑上面的例子:

首先掷一次骰子,结果为1,此时对应的最大概率骰子序列就是D4,因为D4产生1的概率是\(\frac{1}{4}\),高于\(\frac{1}{6}\)和\(\frac{1}{8}\)。

把这个情况扩展,掷两次骰子,结果分别为1,6。这是问题变得复杂起来,要计算3个值,分别是第二个骰子是D4,D6,D8的最大概率。显然,要取得最大概率,第一个骰子必须为D4,此时,第二个骰子取到D6的最大概率如下:

\(P2(D6)\) \(=P(D4)*(1|D4)*(D6|D4)*P(6|D6)\) \(=\frac{1}{3}*\frac{1}{4}*\frac{1}{3}*\frac{1}{6}=\frac{1}{216}\)

同样的,可以计算第二个骰子是D4或D8时的最大概率。发现,第二个骰子取到D6的概率最大。而使这个概率最大时,第一个骰子为D4。所以最大概率骰子序列就是D4,D6。

继续扩展,掷三次骰子,计算第三个骰子分别是D4,D6,D8的最大概率。再次发现,要取到最大概率,第二个骰子必须为D6。这时,第三个骰子取到D4的最大概率如下:

\(P3(D4)\) \(=P2(D6)*P(D4|D6)*P(3|D4)\) \(=\frac{1}{216}*\frac{1}{3}*\frac{1}{4}=\frac{1}{2592}\)

和计算两个骰子序列概率的方法一样,还可以计算第三个骰子是D6或D8时的最大概率。可以发现,第三个骰子取D4的概率最大。而使这个概率最大时,第二个骰子为D6,第一个骰子为D4。

既然掷骰子1到3次可以算,掷多少次都可以,以此类推。

  • 可以发现,要求最大概率骰子序列时要做下面几件事情:

首先,不管序列多长,要从序列长度为1算起,算序列的长度为1时取到每个骰子的最大概率。然后,逐渐增加长度,每增加一次长度,重新算一遍在这个长度下最后一个位置取到每个骰子的最大概率。因为上一个长度下取到每个骰子的最大概率都算过了,重新计算其实不难。当算到最后一位时,就知道最后一位哪个骰子的概率最大了。然后,要把对应这个最大概率的序列从后往前推出来。这就是在刚刚掷骰子的例子中展示出的完整维比特算法。


维特比算法应用——打字提示功能

维比特算法研究的是一种链的可能性问题。现在应用最广的领域是打字提示功能。

Windows中能够使用的输入法有很多种,在这里以全拼为例。在使用输入法时,输入的是英文字母,在这个应用中,隐藏的序列时真实要输入的中文字符和词汇,显示的部分是输入的英文字符。

在输入jin时,输入法软件会猜测想要输入“近”,“斤”,“今”等。而这种排序不是瞎猜的,通常是根据统计数据而来。也就是说“近”,“斤”,“今”这样的顺序一般是根据使用人的输入习惯形成的——在平时打字聊天的过程中“近”出现在词汇或句子输入开始的概率大于“斤”,而“斤”大于“今。

但是输入了tian时就不一样了。“今天”作为一个词汇,比其他任何一个被拼作jintian的词汇都使用得更为高频。也可以理解为,当输入tian时,由jin(今)到tian(天)这条路径的概率是最高的,这是把“今天”这个词汇放在第一个的原因。

后面输入了几个其他的完整词汇:“我们”,“应该”,“做些”,“什么”,输入法也会继续对这些词汇在句子中的路径概率进行计算,每次输入都会猜测一次到目前的输入状态为止最有可能的那条路径,那么看到的这个第一顺位的词汇,准确说是一个句子——“今天我们应该做些什么”就是猜测到的最优的结果,它比任何一种路径产生的概率都要高。

在具体操作的时候,这个马尔可夫模型的训练(拼音串的输入与最终产生汉字串的输出)应该来源于本地的输入者的习惯和互联网的结合。

下面模拟输入法的猜测方法给出一个算法示例,先给出各级转换矩阵,如下表:

  • 表1 jin概率矩阵:

[mdx_table header=”true”]
jin
概率
—–

0.3
—–

0.2
—–

0.1
—–

0.06
—–

0.03
[/mdx_table]

  • 表2 jin-tian转移矩阵:

[mdx_table header=”true”]
jin-tian





—–

0.001
0.001
0.001
0.001
0.001
—–

0.001
0.001
0.001
0.001
0.001
—–

0.990
0.001
0.001
0.001
0.001
—–

0.001
0.001
0.850
0.001
0.001
—–

0.001
0.001
0.001
0.001
0.001
[/mdx_table]

  • 表3 wo概率矩阵:

[mdx_table header=”true”]
wo
概率
—–

0.400
—–

0.150
—–

0.090
—–

0.050
—–

0.030
[/mdx_table]

  • 表4 wo-men转移矩阵:

[mdx_table header=”true”]
wo—men





—–

0.970
0.001
0.003
0.001
0.001
—–

0.001
0.001
0.001
0.001
0.001
—–

0.001
0.001
0.001
0.001
0.001
—–

0.001
0.001
0.001
0.001
0.001
—–

0.001
0.001
0.001
0.001
0.001
[/mdx_table]

在输入一个完整拼音后,用户就会按回车或数字把输入备选框中的汉字输出,当单字词输出时就会由统计产生“jin概率矩阵”和“wo概率矩阵”这样的统计结果,而当用户输入的是一个双字词是就会产生“jin-tian转移矩阵”和“wo-men转移矩阵”这样的统计结果。而且每个双字词,三字词等的输入统计都用这种方法。在输入双字词汉字拼音时会根据转移概率表进行计算。多个词相连就是多个转移矩阵的概率相乘计算,从而得到概率最大的输入可能项。

示例代码如下:

#维特比算法模拟输入法
import numpy as np

jin=['近','斤','今','金','尽']
jin_per=[0.3,0.2,0.1,0.06,0.03]

jintian=['天','填','田','甜','添']
jintian_per=[
[0.001,0.001,0.001,0.001,0.001],
[0.001,0.001,0.001,0.001,0.001],
[0.099,0.001,0.001,0.001,0.001],
[0.002,0.001,0.085,0.001,0.001],
[0.001,0.001,0.001,0.001,0.001]
]

wo=['我','窝','喔','握','卧']
wo_per=[0.400,0.150,0.090,0.050,0.030]

women=['们','门','闷','焖','扪']
women_per=[
[0.970,0.001,0.003,0.001,0.001],
[0.001,0.001,0.001,0.001,0.001],
[0.001,0.001,0.001,0.001,0.001],
[0.001,0.001,0.001,0.001,0.001],
[0.001,0.001,0.001,0.001,0.001]
]

N=5

def found_from_oneword(oneword_per):
    index=[]
    values=[]
    a=np.array(oneword_per)
    for v in np.argsort(a)[::-1][:N]:
        index.append(v)
        values.append(oneword_per[v])
    return index,values

def found_from_twoword(oneword_per,twoword_per):
    last=0
    for i in range(len(oneword_per)):
        current=np.multiply(oneword_per[i],twoword_per[i])
        if i==0:
            last=current
        else:
            last=np.concatenate((last,current),axis=0)
    index=[]
    values=[]
    for v in np.argsort(last)[::-1][:N]:
        index.append([v//5,v%5])
        values.append(last[v])
    return index,values

def perdict(word):
    if word=='jin':
        for i in found_from_oneword(jin_per)[0]:
            print(jin[i])
    elif word=='jintian':
        for i in found_from_twoword(jin_per,jintian_per)[0]:
            print(jin[i[0]]+jintian[i[1]])
    elif word=='wo':
        for i in found_from_oneword(wo_per)[0]:
            print(wo[i])
    elif word=='women':
        for i in found_from_twoword(wo_per,women_per)[0]:
            print(wo[i[0]]+women[i[1]])
    elif word=='jintianwo':
        index1,values1=found_from_oneword(wo_per)
        index2,values2=found_from_twoword(jin_per,jintian_per)
        last=np.multiply(values1,values2)
        for i in np.argsort(last)[::-1][:N]:
            print(jin[index2[i][0]]+jintian[index2[i][1]]+wo[i])
    elif word=='jintianwomen':
        index1,values1=found_from_twoword(jin_per,jintian_per)
        index2,values2=found_from_twoword(wo_per,women_per)
        last=np.multiply(values1,values2)
        for i in np.argsort(last)[::-1][:N]:
            print(jin[index1[i][0]]+jintian[index1[i][1]]+wo[index2[i][0]]+women[index2[i][1]])
    else:
        pass

if __name__=='__main__':
    print('jin:')
    perdict('jin')
    print('jintian:')
    perdict('jintian')
    print('wo:')
    perdict('wo')
    print('women:')
    perdict('women')
    print('jintianwo:')
    perdict('jintianwo')
    print('jintianwomen:')
    perdict('jintianwomen')

根据算法,示例输入“jin”,“jintian”,“wo”,“women”,“jintianwo”和“jintianwomen”时的排序为:
[mdx_fold title=”输出结果:”]
jin:





jintian:
今天
金田
近天
近填
近田
wo:





women:
我们
我闷
我门
我焖
我扪
jintianwo:
今天我
金田窝
近天喔
近填握
近田卧
jintianwomen:
今天我们
金田我闷
近田我扪
近填我焖
近天我门
[/mdx_fold]

在实际应用过程中,这个概率矩阵会是一个系数矩阵。而且转移概率足够小,如小于0.001时可以认为是输入统计中的噪声点,不进行词汇输入推荐,这样输入备选框的前面也只会出现高频输入词汇,这样备选框比较简洁。


想了解更多关于大数据和机器学习:

分类
程序和算法

大数据与机器学习 基础篇 分类 决策树

分类算法是机器学习中的一个重点,也是人们常说的“有监督的学习”。这是一种利用一系列已知类别的样本来对模型进行训练调整分类器的参数,使其达到所要求性能的过程,也成为监督训练或有教师学习。

注:本文中用到的Python及其模块安装教程参见


决策树

决策树是分类时第一种常用的方式,这种方式几乎可以无师自通。

下面举一个例子:

假如某大龄女青年在相亲网站上进行海选,因为资源太多而自己精力有限,所以肯定是要进行相亲决策的,如下图所示:

  • 相亲决策树:

可以看出,这个人的决策策略为:

  • 年龄35以上的直接拉黑,年龄35以下可以考录见面。
  • 年收入20万元以上的属于比较有能力的高质量男性,其他条件可以放宽。
  • 如果学历为硕士以上,身高不够175cm也可以;如果学历是硕士及以下,那么身高必须在175cm以上。
  • 如果年收入在20万元以下,考虑其他因素。
  • 如果学历为硕士以上,身高不够180cm也可以;如果学历是硕士及以下,那么身高必须在180cm以上。

本例是根据样本——对该大龄女青年打招呼的男性的情况从树根开始走决策树,最后决定是相亲还是不相亲。

这种决策树的归纳过程是一个“自底而上”的认识过程。

例如她的相亲策略不是由她直接口述得来的,而是她在网站上不断地相亲,通过计算机学习和归纳总结出来的。

况且更多时候我们是没有机会听别人陈述的,更多的是在观察和总结中认识世界,而且这种归纳总结出来的过程可能比陈述的过程更为客观和准确——所谓察其言不如观其行。

下面继续分析上述相亲的例子,由她之前的相亲信息表归纳出她的相亲策略。

  • 相亲信息表:

[mdx_table header=”true”]
网站ID
年龄(岁)
身高(cm)
年收入(万元)
学历
是否相亲
—–
XXXXXX
25
179
15
大专
N
—–
XXXXXX
33
190
19
大专
Y
—–
XXXXXX
28
180
18
硕士
Y
—–
XXXXXX
25
178
18
硕士
Y
—–
XXXXXX
46
177
100
硕士
N
—–
XXXXXX
40
170
70
本科
N
—–
XXXXXX
34
174
20
硕士
Y
—–
XXXXXX
36
181
55
本科
N
—–
XXXXXX
35
170
25
硕士
Y
—–
XXXXXX
30
180
35
本科
Y
—–
XXXXXX
28
174
30
本科
N
—–
XXXXXX
29
176
36
本科
Y
[/mdx_table]

这里我们考虑一个问题,为什么上图中以年龄与35岁相比作为树根。试想一下,其他的数据项能不能做树根?另外,为什么要以35岁作为树根分裂的条件呢,为什么不是34岁或36岁?对于是否存在一种比较科学或客观的方法能够找到这个描述最简洁的方式,下面引出一个新概念——信息增益。

信息增益

首先这里引用信息熵的概念。

信息熵是用来描述信息混乱程度或者说确定程度的一个值。混乱程度越高,熵越大;混乱程度越低,熵越小。

整个样本集合的熵如下:

\(Info=-\sum_{i=1}^mp_ilog_2p_i\)

关于信息熵的具体解释参见:大数据与机器学习 入门篇

在这里,m的数量就是最后分类(决策)的种类,在本例中m=2——要么见面要么不见面。

\(p_i\)是这个决策项产生的概率。本例中有两个决策项,一个是N(不相亲),概率为\(\frac{5}{12}\),一个是Y(相亲),概率为\(\frac{7}{12}\),这个概率就是从拿到的完整的相亲记录里得到的结论,信息熵为:

\(-\bigl(\frac{5}{12}log_2\frac{5}{12}+\frac{7}{12}log_2\frac{7}{12}\bigr)=0.98bit\)

这里的信息熵也叫期望信息。

期望信息的含义就是:我们要得到这些信息,才能消除这个问题的不确定性。即要做出最后决策需要0.98bit的信息。从熵的定义也不难看出,熵越大说明信息的混乱程度越高,需要更多的信息才能做好决策。

而我们现在要挑选的“树根”,挑选的原则就是要尽可能地消除不确定性,最好做出“一刀切”,即一刀下去就把两个类分清楚。

接下来要考虑的事情就是,以哪个字段做树根能够使得消除信息混杂的能力最强。

假设用某一字段A来划分,在这种划分规则下的信息熵为:

\(Info_A=-\sum_{j=1}^vp_j·Info(A_j)\)

例如,取A为学历,这里v=3,表示有3组,\(A_j\)分别表示大专,本科,硕士,\(p_j\)为各种分组产生的概率,\(Info(A_j)\)是当前分组状态下的期望信息值。具体计算学历的信息熵为:

\(Info_A\) \(=-[(大专项)+(本科项)+(硕士项)]\) \(=-[(p_{大专}×大专分隔熵)+(p_{本科}×本科分隔熵)+(p_{硕士}×硕士分隔熵)]\) \(=-[\frac{2}{12}×(\frac{1}{2}·log_2\frac{1}{2}+\frac{1}{2}·log_2\frac{1}{2})]-[\frac{5}{12}×(\frac{3}{5}·log_2\frac{3}{5}+\frac{2}{5}·log_2\frac{2}{5})]-[\frac{5}{12}×(\frac{4}{5}·log_2\frac{4}{5}+\frac{1}{5}·log_2\frac{1}{5})]\) \(≈0.872\)

信息增益如下:

\(Gain(学历)\)
\(=Info-Info_{学历}\)
\(=0.98-0.872\)
\(=0.108bit\)

这就是用“学历”为根的信息增益,即知道了学历信息后所能消除的不确定性大小。

枚举类型字段的期望信息计算方法:

#枚举类型字段的期望信息的计算方法
import math

#学历分类中大专,本科,硕士占比
education=(2.0/12,5.0/12,5.0/12)
#大专分类中相亲占比
junior_college=(1.0/2,1.0/2)
#本科分类中相亲占比
undergraduate=(3.0/5,2.0/5)
#硕士分类中相亲占比
master=(4.0/5,1.0/5)
#学历各分类中相亲占比
date_per=(junior_college,undergraduate,master)

#“相亲”字段划分规则下的熵
def info_date(p):
    info=0
    for v in p:
        info+=v*math.log(v,2)
    return info

#使用“学历”字段划分规则下的熵
def infoA():
    info=0
    for i in range(len(education)):
        info+=-education[i]*info_date(date_per[i])
    return info

print(infoA())

在归纳决策树的时候,要计算所有字段的信息增益,选取信息增益最大的字段作为决策树根。

连续型变量

计算一下用“年龄”字段做树根能否得到最大的信息增益。

但是“年龄”字段比较麻烦,是一个连续型的字段,不像学历,只有3个枚举值。这种时候通常是在这个字段上找到一个最佳分裂点,然后一刀切下去,让它的信息增益最大。

例如字段中的年龄值,做一个排序:

(25,25,28,28,29,30,33,34,35,36,40,46)

我们选择相邻值的中间值作为分裂点,从中位点开始,向两边计算,每次切割产生一个信息熵,当从m点分裂的信息熵比从m-1和m+1点的信息熵都小时,我们就认为找到了这个点。

连续类型字段的期望信息计算方法:

#连续类型字段的期望信息计算方法
import numpy as np

#年龄
age=[25,25,28,28,29,30,33,34,35,36,40,46]
#是否相亲:1:是,0:否
date=[0,1,1,0,1,1,1,1,1,0,0,0]

splits=[]
for i in range(1,len(age)):
    splits.append((age[i]+age[i-1])/2)

infoAs=[]

age=np.array(age)
date=np.array(date)

#“相亲”字段划分规则下的熵
def info_date(p):
    info=0
    for v in p:
        if v==0:
            return 0
        info+=v*np.math.log(v,2)
    return info

def infoA():
    info=0
    for i in range(len(split_per)):
        info+=-split_per[i]*info_date(date_per[i])
    return info

for i in range(len(splits)):
    #以split为分裂点
    split=splits[i]

    left=len(age[age<=split])

    #左,右分类中的数量占总数的百分比
    split_per=(left/len(age),1-left/len(age))

    #左边分类中相亲占比
    date_left=(len(date[:left][date[:left]==1])/len(date[:left]),1-len(date[:left][date[:left]==1])/len(date[:left]))

    #右边分类中相亲占比
    date_right=(len(date[left:][date[left:]==1])/len(date[left:]),1-len(date[left:][date[left:]==1])/len(date[left:]))

    #左,右各分类中相亲占比
    date_per=(date_left,date_right)

    infoAs.append(infoA())

for i in range(len(splits)):
    print(splits[i],' : ',infoAs[i])

[mdx_fold title=”输出结果为:”]
25.0 : 0.9757921620455572
26.5 : 0.9757921620455572
28.0 : 0.9696226686166434
28.5 : 0.9696226686166434
29.5 : 0.979279160376092
31.5 : 0.9591479170272448
33.5 : 0.90804974601998
34.5 : 0.8112781244591329
35.5 : 0.5731533798814652
38.0 : 0.7344090826922439
43.0 : 0.8668552792172535
[/mdx_fold]

由输出结果可得,选取35.5作为分裂点信息增益最大。

总结一下

最后归纳总结一下构造整棵树时的思路,应该是遵循下面这样的方式:

  • 第一步,找到信息增益最大的字段A和信息增益最大的切分点v(不管是连续类型还是枚举类型)。
  • 第二步,决定根节点的字段A和切分点v。
  • 第三步,把字段A从所有的待选字段列表中拿走,再从第一步开始找。这时决策虽然已经走了一步了,但是根节点也分裂成了两个分支,那么每一个分支各自又形成一个完整的决策树的选择过程。不同的是:可选的字段不一样了,因为A字段被去掉了;此外,在这个分支上的样本也比原来少了,因为两个分支分隔了整个样本,使得一个部分分支只拥有样本的一部分。

减枝法

在决策树的构建过程中,还可以使用“减枝法”进行树的修剪,分为“前减枝”和“后减枝”两种方法。

  • 前减枝:提前终止树的构造,如只用了两个字段,两层树就已经构造出完整的整个树了,那么就可以提前终止树的构造,保持了树的精简性。
  • 后减枝:在决策树完全构造完之后,如建模一共使用了7个字段,全部用上,这样就形成了一个7层的树,如果一个分支下分类已经比较“纯粹”了,就没必要在通过其他条件分支来进行细化,那么整个枝可以直接减掉变成一个叶。

随机森林

随机森林算法是一种并行性比较好的算法规则。

与决策树归纳的过程类似,随机森林也是一个构造决策树的过程,只不过它不是要构造一棵树,而是构造许多棵树。

在决策树的构造过程中会遇到过拟和欠拟的问题,但是在随机森林算法中,通常在一棵树上是不会追求及其精确的拟合的,相反,希望的是决策树的简洁和计算的快速。

随机森林算法的步骤和原则如下:

  • 第一步,随机挑选一个字段构造树的第一层。
  • 第二步,随机挑选一个字段构造树的第二层。
  • ……
  • 第n步,随机挑选一个字段构造树的第n层。
  • 在本棵树构造完毕之后,还需要照相同的方式建造m棵决策树。

补充原则如下:

  • 每棵树的层级通常都比较浅。
  • 每棵树的分类都不能保证分类精度很高。
  • 一个样本进行分类时同时对这m棵决策树做分类概率判断。

人们会为一个训练集构造若干棵决策树,通常可能是几十甚至上百棵,具体会根据样本属性的数量和杂乱程度来决定。当有新样本需要进行分类时,同时把这个样本给这几棵树,然后用民主投票的方式来决定新样本应该属于哪个类,哪一类得票多就归为哪一类。

下面把上述例子用随机森林的方式来实现:

#随机森林
from sklearn.ensemble import RandomForestClassifier

#学历:0:大专,1:本科,2:硕士
#年龄,身高,年收入,学历
X=[
[25,179,15,0],
[33,190,19,0],
[28,180,18,2],
[25,178,18,2],
[46,177,100,2],
[40,170,70,1],
[34,174,20,2],
[36,181,55,1],
[35,170,25,2],
[30,180,35,1],
[28,174,30,1],
[29,176,36,1]
]

#是否相亲:0:N,1:Y
y=[0,1,1,1,0,0,1,0,1,1,0,1]

#现在把训练数据和对应的分类放入分类器中进行训练
clf=RandomForestClassifier().fit(X,y)

#预测下面的人是否相亲
p=[[30,185,18,2],[30,170,70,0]]

print(clf.predict(p))

预测接下来的两个人是否安排相亲,分别为:

年龄30,身高185,年收入18万元,硕士学历,结果为安排相亲;
年龄30,身高170,年收入70万元,大专学历,结果为不安排相亲。


想了解更多关于大数据和机器学习:

分类
程序和算法

C++实现哈希表词频统计

还是上学期数据结构课的练习,这里的代码写得还不是很好,不过可以作为理解哈希表的参考(๑•̀ㅂ•́)و✧

实现功能

1.文件输入输出

1)将文件按单词输入
2)将单词及统计结果保存到文件

2.单词的处理

1)将大小写不同的相同单词合并
2)删除多余的标点符号

3.单词词频统计

算法简述

算法流程图

1.流程图说明说明:

1)算法中采用一个结构体存储单词和词频,结构体包含一个string用于存储单词,一个int用于存储词频,以及一个指向下一节点的指针用于构建链表。

struct hashNode
{
    string word;
    int count;
    hashNode *next;
};

2)哈希地址表是一个数组,下标范围覆盖了哈希函数的整个值域。数组中存储的是指向链表节点(存储单词和词频的结构体)的指针。地址表被声明时初始化为空。

hashNode *index[range] = { NULL };

3)所有单词统计完成后哈希表即失去作用,从链表把数据转入vector后输出。

2.流程图

各部分算法

1.字符串处理

从文本中读取的字符串有可能包含空格以及标点符号,还会存在首字母大写的问题。空格已经在输入输出中由ifstream自动忽略,而标点符号和大小写问题则可以通过ASCII码范围过滤解决。

string handleString(string inputWord)
{
    string word;
    for (int i = 0; i < inputWord.size(); i++)
    {
        if (inputWord[i] >= 65 && inputWord[i] <= 90)
        {
            char tmp = inputWord[i] + 32;
            word.push_back(tmp);
        }
        else if (inputWord[i] >= 97 && inputWord[i] <= 122)
        {
            word.push_back(inputWord[i]);
        }
        else if (inputWord[i] == 45 || inputWord[i] == 39)
        {
            word.push_back(inputWord[i]);
        }
    }
    return word;
}

2.计算哈希值

为了使哈希表尽可能达到最优性能,要使得碰撞概率尽可能小,也就是哈希函数值域和随机性足够大;但是又因为分配过大的内存也会对性能造成影响,故要选择合理的范围。

在本程序实现中,我考虑使用0到9972作为哈希函数值域范围,因为9973是10000内最大的素数,在取模运算中同模概率低;同时它也足够小,因为一个包含9973个元素的指针数组不至于影响计算机性能。

算法上,本实现中使用字符串(单词)各个字母的ASCII码进行迭代运算。ASCII码是将char与int建立对应关系的最简便映射,虽然其随机性远不如SHA1和MD5等高强度哈希算法所采用的对二进制进行操作,但毕竟词频统计的实际应用场景中ASCII码足以保证足够低的碰撞率。

每一位(每一个字母ASCII码)参与计算后的结果将继续参与运算。这种思路是参考了大多数伪随机数算法,通常以系统时间等变量取得第一个随机种子后,之后的随机种子会用第一个随机数运算结果。在计算哈希值中采用这种方法,理论上可以提高哈希值离散度。

int Hash::calculateHash(string word)
{
    int h = 0;
    for (int i = 0; i < word.size(); i++)
    {
        h = 31 * h + word[i];
        h %= 65536;
    }
    int key = h%range;
    return key;
}

在循环中会对65536取模,这是考虑到当单词较长时可能发生计算结果过大影响性能甚至超出int能表示的范围的问题。通过这步取模,可以保证结果不太大,不影响性能。

最后对数组大小9973(已定义为range)取模,保证哈希值值域在数组下标范围内。

3.碰撞检测和计数

当哈希值对应的哈希地址表地址为空时,就表示该单词此前未出现,可以直接添加链表节点开始统计。当该地址非空时,要先判断是发生了碰撞还是找到单词,如果找到单词则计数自增,如果发生碰撞则向后遍历各节点并循环进行判断直到找到单词或地址为空。这种遍历避免冲撞的方法是最简单也是执行效率最低的,但是考虑到哈希函数冲撞率本身很低,发生的冲撞不至于影响整体性能,所以在这一实现中我才用最简单的遍历算法解决冲撞。

void Hash::hashWord(string word)
{
    int key = calculateHash(word);
    hashNode *tmp = index[key];
    while (tmp != NULL)
    {
        if (tmp->word == word)
        {
            tmp->count++;
            return;
        }
        tmp = tmp->next;
        collision += 1;
    }
    hashNode *newNode = new hashNode;
    newNode->count = 1;
    newNode->word = word;
    newNode->next = index[key];
    index[key] = newNode;
}

当发生冲撞时,对冲撞计数的变量自增,便于运行结束时输出。

代码

1.main.cpp

#include <string>
#include <fstream>
#include <iostream>
#include <iomanip>
#include "Hash.h"
#include "Globle.h"

using namespace std;

int main()
{
    string input, inputFileName, outputFileName;
    cout << "请输入运算命令的完整文件名(包含路径和扩展名):";
    cin >> inputFileName;
    cout << "请输入希望保存运算结果的文件名(包含路径和扩展名):";
    cin >> outputFileName;
    ifstream inputFile;
    ofstream outputFile;
    inputFile.open(inputFileName);
    outputFile.open(outputFileName);
    string inputWord;
    Hash hashTable;
    clock_t start, end;
    double time;
    start = clock();
    while (inputFile >> inputWord)
    {
        hashTable.hashWord(handleString(inputWord));
    }
    inputFile.close();
    hashTable.writeOutTable();
    end = clock();
    for (int i = 0; i < hashTable.outTable.size(); i++)
    {
        outputFile << setw(30) << setiosflags(ios::left) << hashTable.outTable[i].word << ' ' << hashTable.outTable[i].count << endl;
    }
    time = (double)(end - start);
    outputFile << "程序主体部分执行时间(毫秒):" << time;
    outputFile.close();
    return 0;
}

2.Hash.h

#pragma once
#include <vector>
#include <string>
#include "Globle.h"

class Hash
{
    hashNode *index[range] = { NULL }; 
    int collision = 0;
public:
    Hash();
    ~Hash();
    void hashWord(string word);
    void writeOutTable();
    int calculateHash(string word);
    vector <hashNode> outTable;
};

3.Hash.cpp

#include "Hash.h"
#include "Globle.h"
#include <string>

Hash::Hash()
{
}

Hash::~Hash()
{
}

void Hash::hashWord(string word)
{
    int key = calculateHash(word);
    hashNode *tmp = index[key];
    while (tmp != NULL)
    {
        if (tmp->word == word)
        {
            tmp->count++;
            return;
        }
        tmp = tmp->next;
        collision += 1;
    }
    hashNode *newNode = new hashNode;
    newNode->count = 1;
    newNode->word = word;
    newNode->next = index[key];
    index[key] = newNode;
}

int Hash::calculateHash(string word)
{
    int h = 0;
    for (int i = 0; i < word.size(); i++)
    {
        h = 31 * h + word[i];
        h %= 65536;
    }
    int key = h%range;
    return key;
}

void Hash::writeOutTable()
{
    for (int i = 0; i < range; i++)
    {
        for (hashNode *tmp = index[i]; tmp != NULL; tmp = tmp->next)
        {
            hashNode outTemp;
            outTemp.count = tmp->count;
            outTemp.word = tmp->word;
            outTemp.next = NULL;
            outTable.push_back(outTemp);
        }
    }
}

4.Globle.h

#pragma once
#include <string>
#include <vector>

#define range 9973//¹þÏ£±í´óС
#define multiply 23//¶¨Òå³Ë·¨ÒòÊý

using namespace std;

string handleString(string inputWord);

struct hashNode
{
    string word;
    int count;
    hashNode *next;
};

5.Globle.cpp

#include "Globle.h"

string handleString(string inputWord)
{
    string word;
    for (int i = 0; i < inputWord.size(); i++)
    {
        if (inputWord[i] >= 65 && inputWord[i] <= 90)
        {
            char tmp = inputWord[i] + 32;
            word.push_back(tmp);
        }
        else if (inputWord[i] >= 97 && inputWord[i] <= 122)
        {
            word.push_back(inputWord[i]);
        }
        else if (inputWord[i] == 45 || inputWord[i] == 39)
        {
            word.push_back(inputWord[i]);
        }
    }
    return word;
}
分类
程序和算法

大数据与机器学习 基础篇 分类 朴素贝叶斯

分类算法是机器学习中的一个重点,也是人们常说的“有监督的学习”。这是一种利用一系列已知类别的样本来对模型进行训练调整分类器的参数,使其达到所要求性能的过程,也成为监督训练或有教师学习。

注:本文中用到的Python及其模块安装教程参见


朴素贝叶斯

贝叶斯决策理论是主观贝叶斯派归纳理论的重要组成部分。

贝叶斯决策理论方法是统计模型决策中的一个基本方法,基本思想如下:

  1. 已知类条件概率参数表达式和先验概率。
  2. 利用贝叶斯公式转换成后验概率。
  3. 根据后验概率大小进行决策分析。

我们先给出贝叶斯公式:

  • 设\(D_1,D_2,……,D_n\)为样本空间S的一个划分,如果以\(P(D_i)\)表示\(D_i\)发生的概率,且\(P(D_i)>0(i=1,2,……,n)\)。对于任何一个事件\(x\),\(P(x)>0\),则有:
\( P(D_j|x)=\frac{P(x|D_j)P(D_j)}{\sum_{i=1}^nP(x|D_i)P(D_i)}\)

最后等式两边都化简为:

\( P(D_j|x)P(x)=P(x|D_j)P(D_j)\)

贝叶斯分类器通常是基于这样一个假定:“给定目标值的属性之间相互条件独立”。

基于这种“朴素”的假定,贝叶斯公式一般简写为:

\(P(A|B)P(B)=P(B|A)P(A)\)

也称为朴素贝叶斯公式——Naive Bayesian。

  • \(P(A)\)叫做\(A\)事件的先验概率,就是一般情况下,认为\(A\)发生的概率。
  • \(P(B|A)\)叫做似然度,是假设\(A\)条件成立的情况下发生\(B\)的概率。
  • \(P(A|B)\)叫做后验概率,在\(B\)发生的情况下发生\(A\)的概率,也就是要计算的概率。
  • \(P(B)\)叫做标准化常量,和\(A\)的先验概率定义类似,就是一般情况下,\(B\)的发生概率。

天气的预测

天气预报一般是基于天气图,卫星云图和雷达图来做的,预报的内容很多,如降水情况,风力,风向等。

大气变化是混沌的(Chaotic)——这一瞬间的天气情况会由于之后一系列不可预测的扰动行为而变得不可预测。也可以说,即便两个完全一模一样的瞬时大气状况,由于之后一系列不可预测的行为不断进行干扰而产生完全两种截然不同的天气结果。

在长期的研究天气变化的过程中,会发现极端的变化还是少数,大部分的天气变化有一定的规律可循。虽然确实不知道未来大气变化的具体情况,但是还是有一些经验可以借鉴的。如今天天气非常阴沉,那下雨的可能性就是比晴天的时候大很多。

利用统计学知识和手段,从大量的历史数据中可以尝试做出一个贝叶斯分类模型,从而判断出降水概率——依据现在的大气状况,未来24小时降水概率为多少。

为了简化说明过程,示意性只列出一个地区10天中每一天的天气状况。

  • 一个地区10天中每天的天气情况表:

[mdx_table header=”true”]
日期
天气情况
—–
1

—–
2

—–
3
降水
—–
4

—–
5
降水
—–
6

—–
7

—–
8
多云
—–
9
降水
—–
10
降水
[/mdx_table]

根据这10天的天气情况来看,10天里面降水4次,也就是降水概率是20%。但是如果当天是阴天呢?阴天的日期为2,4两天,它们的第二天都是降水,所以当天是阴天的情况下,降水概率是100%。如果当天是晴天呢?晴天的日期是1,6,7,而晴天的第二天的天气情况分别是阴,晴,多云,降水概率是0%。

在Python的Scikit-learn库中虽然对朴素贝叶斯分类算法做了实现,但是对于建模针对性的问题,分别做了以下几种贝叶斯分类的变种模型封装:

  1. 高斯朴素贝叶斯(Gaussian Naive Bayes);
  2. 多项式朴素贝叶斯(Multinomial Naive Bayes);
  3. 伯努利朴素贝叶斯(Bernoulli Naive Bayes)。

这三种训练的方式非常相近,引用时所写的代码也非常简短。其中,高斯朴素贝叶斯是利用高斯概率密度公式来进行拟合的。多项式朴素贝叶斯多用于高维度向量分类,最常用的场景是文章分类。伯努利朴素贝叶斯一般是这针对布尔类型特征值的向量做分类的过程。

综上所诉,这里采用高斯朴素贝叶斯为宜,下面给出示例代码:

#高斯朴素贝叶斯预测天气
from sklearn.naive_bayes import GaussianNB

#0:晴 1:阴 2:降水 3:多云
weather={0:'晴',1:'阴',2:'降水',3:'多云'}
data_table=[["date","weather"],
[1,0],
[2,1],
[3,2],
[4,1],
[5,2],
[6,0],
[7,0],
[8,3],
[9,2],
[10,2]]

next=len(data_table)

#当天的天气
X=[]
for i in range(1,(len(data_table)-1)):
    X.append([data_table[i][1]])

#当天的天气对应的后一天的天气
y=[]
for i in range(2,len(data_table)):
    y.append(data_table[i][1])

#把训练数据和对应的分类放入分类器中进行训练
for i in range(3):
    clf=GaussianNB().fit(X,y)
    #由最近一天的天气情况预测明天的天气情况
    p=[[y[len(y)-1]]]
    tmp=clf.predict(p)
    tmp=tmp[0]
    X.append([y[len(y)-1]])
    y.append(tmp)
    data_table.append([next,tmp])
    next+=1

for i in range(1,len(data_table)):
    print(data_table[i][0],end=' ')
    print(':',end=' ')
    print(weather[data_table[i][1]])

预测最近三天的天气,输出结果如下:

[mdx_fold title=”预测结果”]
1 : 晴
2 : 阴
3 : 降水
4 : 阴
5 : 降水
6 : 晴
7 : 晴
8 : 多云
9 : 降水
10 : 降水
11 : 降水
12 : 降水
13 : 降水
[/mdx_fold]

由此可见,此模型预测未来3天都将下雨。

为了更好的说明问题,下面我们考虑一个更简单化的模型。

  • 一个地区10天中每天的天气情况表:

[mdx_table header=”true”]
日期
天气情况
—–
1

—–
2

—–
3
降水
—–
4
多云
—–
5

—–
6

—–
7
降水
—–
8
多云
—–
9

—–
10

[/mdx_table]

可以看出这个地区的天气情况具有比较明显的规律性。将天气数据带入上述代码,输出结果为:

[mdx_fold title=”预测结果”]
1 : 晴
2 : 阴
3 : 降水
4 : 多云
5 : 晴
6 : 阴
7 : 降水
8 : 多云
9 : 晴
10 : 阴
11 : 降水
12 : 多云
13 : 晴
14 : 阴
15 : 降水
16 : 多云
17 : 晴
18 : 阴
19 : 降水
20 : 多云
[/mdx_fold]

可以看出,高斯朴素贝叶斯训练很好地找到了样本之间的联系,并准确的预测出了未来10天的天气状况。

在实际的天气预测中,参数要复杂很多,不仅是当天的天气这么一个简单的变量,而且还有温度,湿度,季节,风速等因素,而且统计的日期也是长期的天气样本信息,这些因素的多少直接决定了模型的复杂程度。


疾病的预测

关于罹患疾病的预测——罹患冠心病的概率为20%,罹患阿尔兹海默综合征的概率为5%,罹患帕金森式病的概率为12%等。这些都是通过基因测序得到的预测结论。

这个预测过程也是一个分类过程,训练样本就是大量的个体基因信息和个人的疾病信息。然后通过建模分析,最后得到一个基因片段和罹患疾病之间的概率转换关系,这也是一个比较典型的朴素贝叶斯分类模型。

示意性地做一个训练过程,训练样本如表所示。

  • 训练样本表:

[mdx_table header=”true”]
样本人员编号
基因片段A
基因片段B
高血压
胆结石
—–
1
‘1’
‘1’


—–
2
‘0’
‘0’


—–
3
‘0’
‘1’


—–
4
‘1’
‘0’


—–
5
‘1’
‘1’


—–
6
‘1’
‘0’


—–
7
‘0’
‘1’


—–
8
‘0’
‘0’


—–
9
‘1’
‘0’


—–
10
‘0’
‘1’


[/mdx_table]

[mdx_warning title=”注意”]这种记录是忽略其他建模因素的记录,如这位病患所生活的城市,所从事的工作,是否有烟酒习惯等,这些因素是没有参与建模的,直接假设罹患这些疾病只和基因状况有关。[/mdx_warning]

如果之后有一个用户来做基因测序,测试结果为基因片段A,B分别为1,0,那么他罹患高血压和胆结石这两种疾病的概率分别是多少?

运用朴素贝叶斯公式:

\( P(A|B)=\frac{P(B|A)P(A)}{P(B)}\)
  • 关于公式中的每个值分别代表什么?下面以计算高血压疾病的罹患概率为例:
    1. \(P(A)\)是先验概率,即全局性的高血压罹患概率,在10个人中有3个人患病,即概率为30%。
    2. \(P(B|A)\)是似然度,表示高血压患者中基因是“10”型的患者数量。在3位高血压患者中只有1位患者的基因是“10”型,即概率为33%。
    3. \(P(B)\)是标准化常量,是“10”型基因出现的概率,在10个人中出现了3例,即概率为30%。
    4. \(P(A|B)\)可以代入求解,得到33%。

交叉验证一下,全局中“10”型基因的有3位,即4号,6号和9号。的高血压的只有9号一人,所以高血压的罹患概率是33%。

下面给出训练和分类的代码示例:

#高斯朴素贝叶斯预测罹患概率
from sklearn.naive_bayes import GaussianNB

#基因片段A,基因片段B,高血压,胆结石
#1:是,0:否
data_table=[
[1,1,1,0],
[0,0,0,1],
[0,1,0,0],
[1,0,0,0],
[1,1,0,1],
[1,0,0,1],
[0,1,1,1],
[0,0,0,0],
[1,0,1,0],
[0,1,0,1]]

#基因片段
X=[[1,1],[0,0],[0,1],[1,0],[1,1],[1,0],[0,1],[0,0],[1,0],[0,1]]

#高血压
y1=[1,0,0,0,0,0,1,0,1,0]

#训练
clf=GaussianNB().fit(X,y1)

#预测
p=[[1,0]]
result=[]
result.append(clf.predict(p))

#胆结石
y2=[0,1,0,0,1,1,1,0,0,1]

#巡练
clf=GaussianNB().fit(X,y2)

#预测
result.append(clf.predict(p))

print("基因:",p,"\n罹患高血压:",result[0][0],"\n罹患胆结石:",result[1][0])

输出结果为:

[mdx_fold title=”预测结果”]
基因: [[1, 0]]
罹患高血压: 0
罹患胆结石: 0
即表示基因型为[1,0]的这个用户罹患两种疾病的概率都小于50%。
[/mdx_fold]


小结

在我看来,贝叶斯的理论体系其实揭示的是一种非常典型的人类自身的推测逻辑行为。

例如,在黄昏的时候走在自己居住的小区里,光线很昏暗,前面突然闪过一个影子,从路一边的草丛蹿到另一边,速度较快体型较大,其他的信息都没有铺捉到。这时候大概会猜测,这有可能是一只较大的家犬。而如果是在非洲大草原上,从越野车里同样看到昏暗的草原上蹿过一个速度较快体型较大的动物,也许会猜测那是一头狮子,或一头猎豹。这两种猜测对于铺捉到的对象信息都是非常有限的,而且内容相近,但是得到的却是两种不同的推测。原因很简单,就是因为当时的环境不同,导致两种事件的概率不同,带有比较明确的倾向性。而这种推断的思路就是贝叶斯理论体系的核心内容。


想了解更多关于大数据和机器学习:

分类
程序和算法

大数据与机器学习 基础篇 聚类

聚类(Clustering)指的是一种学习方式,即把物理或抽象对象的集合分组为由彼此类似的对象组成的多个类的分析过程。

注:本文中用到的Python及其模块安装教程参见


K-Means算法

在聚类中K-Means算法是很常用的一个算法,也是基于向量距离来做聚类。

  • 算法步骤如下:
    1. 从n个向量对象中选择任意k个向量作为初始聚类中心。
    2. 根据在步骤1中设置的k个向量(中心对象向量),计算每个对象与这k个中心对象各自的距离。
    3. 对于步骤2中的计算,任何一个向量与这k个向量都有一个距离,有的远有的近,把这个向量的距离它最近的中心向量对象归在一个类簇中。
    4. 重新计算每个类簇的中心对象向量的位置。
    5. 重复3,4两个步骤,直到类簇聚类方案中的向量归类变化极少为止。例如:一次迭代后,只有少于1%的向量还在发生类簇之间的归类漂移,那么就可以认为分类完成。
  • 这里要注意的是:
    1. 需要事先指定类簇的数量。
    2. 需要事先给定初始的类中心。

我们看这个例子:

准备一个中国城市经纬表,一共2259个向量,维度分别为(经度,维度),城市信息已经略去,数据参见附录。

设定类簇数量为5,开始聚类。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans

X=[]
f=open('city.txt')
for v in f:
    X.append([float(v.split(',')[0]),float(v.split(',')[1])])

X=np.array(X)

n_clusters=5

cls=KMeans(n_clusters).fit(X)

cls.labels_

markers=['^','x','o','*','+']
for i in range(n_clusters):
    members=cls.labels_==i
    plt.scatter(X[members,0],X[members,1],s=60,marker=markers[i],c='b',alpha=0.5)

plt.title('China City Area Distribution')
plt.show()
  • 聚类效果如图:

从图片可以清楚地看出,中国的城市被划分成了东北,华中,华南,西部,西北5个区域。

  • K-Means算法所使用的距离计算公式是可以选取的,一般用欧式距离比较简单,也可以用曼哈顿距离。常用的距离度量方法还有余弦相似度。欧式距离度量会受指标不同单位刻度的影响,所以一般需要先进行标准化或者归一化,同时距离越大,个体间差异越大;空间向量余弦夹角的相似度度量不会受指标刻度的影响,余弦值落于区间[-1,1]上,值越大,差异越小。有关余弦相似度的例子在之后的文章中会具体介绍。

层次聚类

层次聚类有两种思路,一种是“凝聚的层次聚类方法”,一种是“分裂的层次聚类方法”。

  • 凝聚的层次聚类方法:在大量的样本中自底而上找那些距离比较近的样本先聚合成小的群,聚合到一定程度再由小的群聚合成更大的群。
  • 分裂的层次聚类方法:先把所有的样本分成若干个大群,再在每个群里各自重新进行聚类划分。

分裂的层次聚类方法可以通过迭代使用K-Means算法来实现,下面重点介绍凝聚的层次聚类方法。

在Scikit-learn库中,提供了一种叫做AgglomerativeClustering的分类算法。首先设置几个观察点散布在整个训练样本中,让这些观察点自下而上不断地进行类簇的合并。这种聚类合并遵循一定的原则,即基于连接度的度量来判断是否要向上继续合并两个类簇。度量有以下3种不同的策略原则:

  1. Ward策略:让所有类簇中的方差最小化。
  2. Maximum策略:力求将类簇之间的距离最大值最小化。
  3. Average linkage策略:力求将类簇之间的距离的平均值最小化。

例如,采用Ward策略:

import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import AgglomerativeClustering

X=[]
f=open('city.txt')
for v in f:
    X.append([float(v.split(',')[0]),float(v.split(',')[1])])

X=np.array(X)

n_clusters=5

cls=AgglomerativeClustering(linkage='ward',n_clusters=n_clusters).fit(X)

cls.labels_

markers=['^','x','o','*','+']
for i in range(n_clusters):
    members=cls.labels_==i
    plt.scatter(X[members,0],X[members,1],s=60,marker=markers[i],c='b',alpha=0.5)

plt.title('China City Area Distribution')
plt.show()
  • 聚类效果如图:


密度聚类

在聚类形状不规则的情况下,K-Means这种基于欧式距离半径的类簇划分方法就不那么好用了,K-Means对于带拐弯,狭长的,不规则形状的类簇的效果就不如对圆形类簇的效果好。

这里使用sklearn中专门用来做密度分类的算法库——sklearn.cluster.DBSCAN。

例如对下表中的国家面积和人口的二维向量进行聚类。

  • 12个国家的面积,人口和GDP报表:

[mdx_table header=”true”]
国家
面积km^2
人口
GDP亿美元
人均GDP美元
—–
中国
9670250
1392358258
99960
7179
—–
印度
2980000
1247923065
18707
1505
—–
美国
9629091
317408015
167997
53101
—–
巴西
8514877
201032714
22429
11311
—–
日本
377873
127270000
49015
38491
—–
澳大利亚
7692024
23540517
15053
64863
—–
加拿大
9984670
34591000
18251
51990
—–
俄罗斯
17098200
143551289
21180
14819
—–
泰国
513115
67041000
3872
5674
—–
柬埔寨
181035
14805358
157
1016
—–
韩国
99600
50400000
12218
24329
—–
朝鲜
120538
24052231
355
1476
[/mdx_table]

import numpy as np
from sklearn.cluster import DBSCAN
import matplotlib.pyplot as plt

X=[]
f=open('country.txt')
for v in f:
    X.append([float(v.split(',')[0]),float(v.split(',')[1])])

X=np.array(X)

#做归一化
a=X[:,:1]/max(X[:,:1])[0]*10000
b=X[:,1:]/max(X[:,1:])[0]*10000
X=np.concatenate((a,b),axis=1)

cls=DBSCAN(eps=2000,min_samples=1).fit(X)

cls.labels_

markers=['^','x','o','*','+']
for i in range(n_clusters):
    members=cls.labels_==i
    plt.scatter(X[members,0],X[members,1],s=60,marker=markers[i],c='b',alpha=0.5)

plt.title("Countrys' Area and population")
plt.show()
  • 聚类效果如图:

  • 需要注意的是:
  1. 归一化:是为了解决由于量纲或单位不同所产生的距离计算问题而进行的权重调整,目的是把两个不同维度的数据都投影到以10000为最大值的正方形区域里。
  2. DBSCAN的参数:
    • eps:设置一个阈值,再根据密度向外扩展的过程中如果发现在这个阈值距离范围内找不到向量,那么就认为这个聚簇已经查找完毕。因为归一化以后的向量都落在10000×10000的区间内,所以这里设置2000作为阈值。
    • min_samples:含义是设置聚簇最小应该拥有多少个向量。例如:如果本题中设置min_samples为3,那么右下角的俄罗斯,上面的中国和印度都将作为噪声点被消除而不会显示。

聚类评估

聚类的质量评估包括以下3个方面:

  1. 估计聚类的趋势。对于给定的数据集,评估该数据集是否存在非随机结构,也就是分布不均匀的情况。如果直接使用各种分类算法套用在样本上,然后返回一些类簇,这些类簇的分布很可能是一种错误的分类,会对人们产生误导。数据中必须存在非随机结构,聚类分析才是有意义的。
  2. 确定数据集中的簇数。上述K-Means算法在一开始就需要确定类簇的数量,并作为参数传递给算法。这里容易让人觉得有点矛盾,即人主观来决定一个类簇的数量的方法是不是可取。
  3. 测量聚类的质量。可以用量化的方法来测量聚类的质量。

聚类趋势

如果样本空间里的样本是随机出现的,本身没有聚类的趋势,那么使用聚类肯定是有问题的。对于聚类的趋势,我们常用霍普金斯统计量(Hopkins Statistics)来进行量化评估。

  1. 从所有的样本向量中随机找到n个向量,把他们成为p向量,每一个向量分别是\(p_1,p_2,……,p_n\)。对每一个向量都在样本空间里找到一个离其最近的向量,然后求距离(欧氏距离),然后用\(x_1,x_2,……,x_n\)来表示这个距离。
  2. 从样本向量的一个子集中随机找到n个向量,把他们成为q向量,每一个向量分别是\(q_1,q_2,……,q_n\)。对每一个向量都在样本空间里找到一个离其最近的向量,然后求距离(欧氏距离),然后用\(y_1,y_2,……,y_n\)来表示这个距离。
  3. 求出霍普金斯统计量H:
\(H=\frac{\sum_{i=1}^ny_i}{\sum_{i=1}^nx_i+\sum_{i=1}^ny_i}\)

如果整个样本空间是一个均匀的,没有聚类趋势(聚簇不明显)的空间,那么H应该为0.5左右。反之,如果是有聚类趋势(聚簇明显)的空间,那么H应该接近于0。

计算霍普金斯统计量:

计算中国城市经纬度的霍布金斯统计量,因为类簇数为5,所以取\(X[np.random.randint(X.shape[0],size=int(X.shape[0]/5))]\)为子集,即随机地选取\(\frac{1}{5}\)个向量作为子集。

代码如下:

import numpy as np

X=[]
f=open('city.txt')
for v in f:
    X.append([float(v.split(',')[0]),float(v.split(',')[1])])

X=np.array(X)

#做归一化
a=X[:,:1]/max(X[:,:1])[0]*100
b=X[:,1:]/max(X[:,1:])[0]*100
X=np.concatenate((a,b),axis=1)

HH=[]

for count in range(10):
    #随机选出10个
    pn=X[np.random.choice(X.shape[0],10,replace=False),:]
    xn=[]
    for i in pn:
        distance_min=1000000
        for j in X:
            if np.array_equal(j,i):
                continue
            distance=np.linalg.norm(j-i)
            if distance_min>distance:
                distance_min=distance
            xn.append(distance_min)

    #在子集X[np.random.randint(X.shape[0],size=int(X.shape[0]/5))]随机选出10个
    XX=X[np.random.randint(X.shape[0],size=int(X.shape[0]/5))]
    qn=XX[np.random.choice(XX.shape[0],10,replace=False),:]
    yn=[]
    for i in qn:
        distance_min=1000000
        for j in XX:
            if np.array_equal(j,i):
                continue
            distance=np.linalg.norm(j-i)
            if distance_min>distance:
                distance_min=distance
            yn.append(distance_min)

    Htmp=float(np.sum(yn))/(np.sum(xn)+np.sum(yn))
    HH.append(Htmp)

print(HH)
H=np.average(HH)
print("H: ",H)

[mdx_fold title=”输出结果为:”]
0.11722419404269789,
0.10871923835253638,
0.09662900489410153,
0.18076604758829162,
0.12209136642523832,
0.16461772084853724,
0.09426555293455162,
0.1924580732471169,
0.14596993096952454,
0.13714731037623953
[/mdx_fold]

平均值为:0.13598884396788355

由此可见,中国城市的经纬度具有明显的聚类趋势。

簇数确定

确定一个样本空间里有多少簇数是很重要的,尤其是在K-Means这种算法里一开始就要求给定要被分成的簇数。况且,簇数的猜测也会影响聚类的结果,簇数太多,样本被分成很多小簇,簇数太少,样本基本没有被分开,都没有意义。

有一种经验的方法就是对于n个样本的空间,设置簇数\(p\)为\(\sqrt{\frac{n}{2}}\),在期望状态下每个簇大约有\(\sqrt{2n}\)个点。但是这种说法只能作为参考。

还有一种方法叫“肘方法”(The Elbow Method),被认为是一种更科学的方式。思路如下:

  • 尝试把样本空间划分为1个类,2个类,……,n个类,要确定哪种分法最为科学,在分成m个类簇的时候会有一个划分方法,在这种划分方法下,每个类簇的内部都有若干个向量,计算这些向量的空间中心点,即计算着m个类簇各自的空间重心在哪里。在计算每个类簇中每个向量和该类簇重心的距离(大于等于0)的和。最后把m个类簇各自的距离和相加得到一个函数var(n),n就是类簇数。

可以想象,这个平方和最大的时刻应该是分1个类——也就是不分类的时候,所有的向量到重心的距离都非常大,这样的距离的和是最大的。那么尝试着划分为2个类,3个类,4个类……随着分类的增多,第m次划分时,每个向量到自己簇的重心的距离会比上一次(m-1)次临近的机会更大,那么这个距离总和就会总体上缩小。极限情况就是最后被分成了n个类簇,n整个空间向量的数量,也就是一个向量一个类簇,每个类簇一个成员。这种情况最后距离的和就变成了0,因为每个向量距离自己(自己就是重心)的距离都是0。

当m从1,2,3……逐步向上增加的过程中,整个曲线的斜率会逐步降低,而且一开始是快速下降的。下降的过程中有一个拐点,这一个点会让人感觉曲线从立陡变成平滑,那么这个点就是要找的点。

这个样本空间被分成m个类簇之后,再分成更多的类簇时,每次的“收获”没有前面每次“收获”那么大,此时的m值就是被认为最合适的类簇数量。这个点在曲线上给人的感觉就像是人的胳膊肘一样,所以被形象地命名为“肘方法”。

例如,对中国城市的经纬表做簇数确定(这里使用曼哈顿距离):

import numpy as np
from sklearn.cluster import KMeans
import matplotlib.pyplot as plt

X=[]
f=open('city.txt')
for v in f:
    X.append([float(v.split(',')[0]),float(v.split(',')[1])])

X=np.array(X)

#做归一化
a=X[:,:1]/max(X[:,:1])[0]*100
b=X[:,1:]/max(X[:,1:])[0]*100
X=np.concatenate((a,b),axis=1)

#曼哈顿距离
def manhattan_distance(x,y):
    return np.sum(abs(x-y))

distance_sum=[]

for N in range(1,10):
    n_clusters=N

    cls=KMeans(n_clusters).fit(X)

    cls.cluster_centers_

    cls.labels_

    distance=0
    for i in range(n_clusters):
        group=cls.labels_==i
        members=X[group,:]
        for v in members:
            distance+=manhattan_distance(np.array(v),cls.cluster_centers_[i])

    distance_sum.append(distance)

for N in range(1,10):
    plt.scatter(N,distance_sum[N-1],s=60,marker='o',c='b',alpha=0.5)

plt.plot(range(1,10),distance_sum,'r',label='Fitted line',)

plt.title("Cluster Number and Distance_sum")
plt.show()
  • 输出结果如图:

这里可以发现选择5个类簇和6个类簇对中国城市经纬表进行聚类都是可以的。

  • 6个类簇的中国城市经纬图(将西北地区又分成了两部分):

有的时候会发现类簇数量可能不是很好确定,这个拐点不清晰。这种时候适当地在计算效率和收益程度上做一个平衡即可,不必太纠结,毕竟聚类是一个无监督的学习。

测定聚类质量

测定聚类质量的方法有很多,一般分为“外在方法”和“内在方法”两种。

所谓外在方法是一种依靠类别基准的方法,即已经有比较严格的类别定义时在讨论聚类是不是足够准确。这里通常使用“BCubed精度”和“BCubed召回率”来进行衡量。但是外在方法适用于有明确的外在类别基准的情况,而聚类是一种无监督的学习,更多的是在不知道基准的状况下进行的,所以我们更倾向于使用“内在方法”。

“内在方法”不会去参考类簇的标准,而是使用轮廓系数(Silhouette Coefficient)进行度量。

计算轮廓系数的思路如下:

  • 对于有n个向量的样本空间,假设它被划分成k个类簇,即\(C_1,C_2,……,C_k\)。对于任何一个样本空间中的向量v来说,可以求一个v到本类簇中其他各点的距离的平均值a(v),还可以求一个v到其他所有各类簇的最小平均距离(即从每个类簇里挑选一个离v最近的向量,然后计算距离),求这些距离的平均值,得到b(v),轮廓系数的定义为:
\(s(v)=\frac{b(v)-a(v)}{max[a(v),b(v)]}\)

一般来说,这个函数的结果在-1和1之间,a(v)表示的是类簇内部的紧凑型,越小越紧凑,而b(v)表示该类簇的其他类簇之间的分离程度。如果函数值接近1,即a(v)比较小而b(v)比较大时,说明包含v的类簇非常紧凑,而且远离其他的类簇。相反,如果函数值为负数,则说明a(v)>b(v),v距离其他类簇比距离自己所在类簇的其他对象更近,那么这种情况就不太好,应该尽可能避免。

为了让聚类中的类簇划分更合理,可以计算类簇中所有对象的轮廓系数的平均值。在一种方案里,如果轮廓系数是负数那么可以直接淘汰,如果是正数则可以在多个方案里进行比较,选择一种轮廓系数接近1的方案。

计算中国城市经纬表聚类后的轮廓系数:

#计算轮廓系数
import numpy as np
from sklearn.cluster import KMeans
import matplotlib.pyplot as plt

X=[]
f=open('city.txt')
for v in f:
    X.append([float(v.split(',')[0]),float(v.split(',')[1])])

X=np.array(X)

#做归一化
a=X[:,:1]/max(X[:,:1])[0]*100
b=X[:,1:]/max(X[:,1:])[0]*100
X=np.concatenate((a,b),axis=1)

#曼哈顿距离
def manhattan_distance(x,y):
    return np.sum(abs(x-y))

n_clusters=6

cls=KMeans(n_clusters).fit(X)

cls.labels_

index=np.random.randint(len(X))

#a(v),X[index]到类簇中其他点的距离的平均值
distance_sum=0
for v in X[cls.labels_==cls.labels_[index]]:
    distance_sum+=manhattan_distance(np.array(X[index]),np.array(v))

av=distance_sum/len(X[cls.labels_==cls.labels_[index]])
print(av)

#b(v),X[index]
distance_sum=0
distance_min=1000000
for i in range(n_clusters):
    if i==cls.labels_[index]:
        continue
    group=cls.labels_==i
    members=X[group]
    for v in members:
        distance=manhattan_distance(np.array(v),X[index])
        if distance_min>distance:
            distance_min=distance
    distance_sum+=distance

bv=distance_sum/(n_clusters-1)
print(bv)

sv=float(bv-av)/max(av,bv)
print(sv)

[mdx_fold title=”输出结果为:”]
5个类簇时的轮廓系数:0.5292298922519052
6个类簇时的轮廓系数:0.6012776282140894
[/mdx_fold]


想了解更多关于大数据和机器学习:

[mdx_fold title=”中国城市经纬表”]
经度,纬度
121.48,31.22
121.24,31.4
121.48,31.41
121.7,31.19
121.76,31.05
121.46,30.92
121.24,31
121.16,30.89
121.1,31.15
121.4,31.73
102.73,25.04
102.48,25.21
102.58,24.68
102.79,24.9
102.44,24.95
103.7,29.32
103.63,28.22
103.91,27.74
104.06,27.61
103.54,27.21
103.97,28.58
104.28,28.08
105.05,27.85
104.86,27.42
102.92,26.9
104.38,28.62
103.79,25.51
104.09,26.24
104.24,25.67
103.97,24.85
103.03,25.35
103.27,26.41
103.82,25.62
104.3,24.88
104.64,25.04
103.12,24.9
103.61,25.41
103.24,24.77
103.25,25.56
102.52,24.35
102.93,24.26
102.75,24.09
102.91,24.68
102.73,24.27
102.15,24.67
102,23.59
101.98,24.06
102.38,24.16
101,22.79
101.03,23.07
100.88,23.9
100.82,24.42
100.71,23.5
101.71,23.4
99.97,22.55
99.47,22.73
101.88,22.58
99.55,22.32
100.09,23.88
100.12,24.44
99.02,23.92
99.25,24.03
99.92,24.58
99.85,23.45
99.24,23.15
99.41,23.56
99.18,25.12
99.15,24.69
98.51,25.01
99.61,24.82
98.7,24.58
100.25,26.86
101.24,26.63
100.76,26.71
100.82,27.29
104.24,23.37
105.09,24.05
104.68,23.42
104.71,23.12
104.4,23.01
104.19,24.03
104.35,23.62
105.6,23.62
102.43,23.35
103.43,24.41
103.41,23.36
102.81,23.17
102.42,23.35
102.48,23.73
103.76,24.52
103.24,22.77
103.23,23.7
102.42,23.01
102.79,23.64
103.98,22.52
103.67,22.68
100.79,22
100.5,21.95
101.56,21.48
101.54,25.01
101.85,25.7
102.36,25.55
102.08,25.15
101.26,25.21
101.34,25.73
101.7,26.07
102.45,25.58
101.58,25.32
101.67,24.68
101.24,25.4
100.24,25.45
99.88,26.53
99.94,26.1
100.55,25.82
100.52,25.34
99.52,25.45
100.18,26.55
100.19,25.69
99.98,25.68
99.39,25.9
100.56,25.48
100.33,25.23
100.51,25.04
98.6,24.41
97.96,24.33
97.93,24.69
98.08,24.08
97.83,24
98.3,24.78
98.82,25.97
98.95,26.55
98.92,26.89
99.29,26.49
98.65,27.73
99.72,27.78
98.93,28.49
99.27,27.15
123.33,25
117.51,15.07
111.65,40.82
111.13,40.72
111.15,40.28
110,40.58
110.52,40.55
110.03,41.03
106.82,39.67
113.08,41.03
113.97,40.88
111.65,39.92
111.42,41.12
112.52,40.93
113.53,41.58
113.15,40.45
112.48,40.52
111.8,40.4
114,41.9
111.68,41.37
111.96,43.65
116.08,43.95
116.48,42.18
114.97,44.03
117.58,44.6
116.97,45.53
113.7,43.85
112.95,42.47
115.3,41.9
115,42.32
116.02,42.25
113.83,42.25
119.73,29.22
117.47,49.58
119.45,49.33
120.08,50.45
121.52,50.8
120.73,49.3
123.5,48.13
122.78,47.98
116.82,48.67
118.23,48.22
123.7,50.58
122.28,43.63
121.32,43.62
121.75,42.72
120.65,42.85
120.87,44.55
118.87,42.28
119.32,41.62
118.02,43.62
118.67,41.95
119.87,42.3
119,42.97
118.65,43.52
119.35,43.98
120.05,43.97
117.48,43.28
110,39.83
110,39.83
111.13,39.68
109.03,38.38
109.77,39.25
107.97,39.12
107.43,38.18
108.7,39.83
110.02,40.42
107.37,40.78
108.28,41.12
106.98,40.33
107.12,40.88
108.52,41.55
108.65,40.75
108.52,40.88
105.68,38.85
101.68,39.2
100.88,41.9
122.08,46.07
121.5,45.4
116.46,39.92
117.1,40.13
116.85,40.37
116.65,40.13
116.67,39.92
116.62,40.32
116.33,39.73
115.98,39.72
125.35,43.88
126.57,43.87
125.15,44.45
125.68,44.52
126.55,44.83
126.83,44.15
125.68,43.53
126.57,43.87
126.97,44.4
127.33,43.75
126.72,42.97
126.03,42.93
129.52,42.93
129.75,43.32
130.35,42.85
129.83,42.98
129,42.52
128.3,42.58
128.18,43.35
125.92,41.49
125.7,40.88
125.65,42.53
126.03,42.68
126.8,42.38
126.4,41.97
127.27,42.33
126.17,41.15
128.17,41.43
124.37,43.17
124.33,43.32
124.82,43.5
125.32,43.33
125.15,42.97
125.5,42.68
123.5,43.52
122.82,45.63
124.18,45.5
124.82,45.2
124.02,45
123.97,44.3
123.13,44.82
122.75,45.35
104.06,30.67
104.32,30.88
104.94,30.57
103.29,30.2
103.86,30.8
104.13,30.82
102.15,26.9
101.56,26.9
103.81,30.97
103.61,31.04
103.94,30.99
104.16,31.13
104.25,30.99
103.78,30.42
103.47,30.42
103.53,30.58
103.69,30.63
104.73,31.48
104.7,31.8
105.21,32.59
104.52,32.42
105.86,32.44
106.33,32.25
105.45,32.03
105.16,31.64
105.06,31.1
105.35,31.23
105.31,30.9
105.58,30.52
105.74,30.78
104.68,31.06
104.37,31.13
104.19,31.32
104.41,31.64
104.44,31.89
105.04,29.59
105.02,30.3
105.3,30.12
104.7,29.57
104.85,29.81
104.6,30.19
104.53,30.38
105.25,29.64
104.56,29.77
104.97,29.24
104.96,28.87
105.06,28.71
105.38,28.77
105.46,28.96
105.78,28.79
105.39,28.91
105.79,28.03
105.44,28.19
104.91,28.6
105.06,28.36
104.81,28.38
104.52,28.4
104.53,28.16
104.15,28.68
103.73,29.59
103.59,29.75
103.38,29.95
103.53,30.04
103.81,29.86
103.81,30.05
103.83,30.22
104.06,29.67
104.09,30
103.93,29.21
103.98,28.96
103.5,29.62
103.53,28.87
103.25,29.23
103.13,29.24
107.36,29.7
107.34,30.36
107.7,29.89
108.13,29.98
108.97,28.47
108.75,28.85
108.81,29.53
108.19,29.29
108.72,29.29
107.13,29.15
108.35,30.83
108.39,31.23
108.67,31.98
109.6,31.42
109.86,31.1
109.52,31.06
108.89,30.99
108.03,30.33
107.78,30.66
106.06,30.8
105.96,31.75
105.97,31.75
106.38,31.52
106.03,31.34
105.84,31.01
106.57,31.07
106.44,31.04
106.61,30.48
106.43,30.55
106.3,30.38
106.74,30.41
107.49,31.23
108.06,32.07
107.71,31.39
107.87,31.1
106.91,30.36
107.21,30.75
106.94,30.85
106.83,32.36
106.73,31.86
107.11,31.59
108.24,31.95
108.18,32
102.97,29.97
102.91,30.17
103.06,30.09
102.81,29.79
102.66,29.4
102.38,29.21
102.78,30.09
102.84,30.36
102.22,31.92
102.55,31.79
101.72,31.93
102.94,33.62
102.95,32.06
103.61,32.64
104.19,33.23
103.61,31.46
103.16,31.42
102.34,30.97
102.03,31.48
100.97,32.3
103.89,31.67
101.95,30.04
100.65,31.38
99.96,31.64
100.28,30.96
98.83,32.23
98.57,31.81
98.06,33.01
100.35,32.3
102.25,29.92
101.87,30.85
101.53,29.01
101,30.03
101.14,30.99
100.28,30.03
99.78,28.93
100.31,29.04
99,30
99.25,28.71
102.29,27.92
102.83,28.03
102.74,28.96
103.62,28.21
102.76,27.07
102.55,26.74
102.21,26.67
102.15,27.4
103.14,28.33
103.22,27.73
102.8,27.7
102.52,27.38
102.42,28.33
102.49,28.66
101.51,27.42
102.15,28.58
101.25,27.9
117.2,39.13
117.83,39.33
116.92,38.93
117.4,40.05
117.3,39.75
117.05,39.4
106.27,38.47
106.24,38.28
106.35,38.55
106.39,39.04
106.54,38.91
106.69,38.82
106.21,37.99
105.94,36.97
106.34,38.1
105.66,37.48
107.41,37.78
105.18,37.51
106.07,38.02
106.28,36.01
105.7,35.97
106.33,35.5
105.64,36.56
106.11,35.63
117.27,31.86
117.16,32.47
116.98,32.62
116.71,32.68
116.77,33.97
116.76,33.92
118.38,31.33
117.82,30.93
117.34,32.93
118.48,31.56
117.03,30.52
116.97,33.63
116.97,33.63
116.34,34.42
116.93,34.19
117.55,33.55
117.89,33.49
117.87,33.14
117.32,33.33
117.19,32.95
118.31,32.33
117.98,32.78
119,32.68
118.44,32.44
118.27,32.1
117.68,32.52
117.4,32.86
117.87,31.62
117.87,31.62
117.47,31.89
118.11,31.7
118.37,31.7
117.75,31.3
117.29,31.23
118.73,31.95
118.49,31.55
119.17,31.14
119.41,30.89
118.41,30.68
118.32,30.91
118.21,31.07
118.95,30.62
117.84,30.64
118.31,29.72
118.19,29.81
118.53,30.28
118.57,30.07
118.44,29.88
117.7,29.86
117.92,29.93
118.13,30.28
117.48,30.19
116.94,31.04
117.21,30.69
116.63,30.41
116.69,30.12
116.13,30.15
116.27,30.42
116.36,30.84
116.53,30.62
116.99,30.08
117.48,30.66
116.49,31.73
116.27,32.32
116.78,32.57
117.15,31.7
116.94,31.45
116.32,31.38
115.87,31.67
115.81,32.89
116.76,33.86
116.21,33.49
116.55,33.25
116.19,33.12
116.26,32.62
115.6,32.63
115.24,33.06
115.34,33.24
115.61,33.16
117,36.65
117.07,36.69
116.73,36.55
117.53,36.72
120.33,36.07
120.42,36.15
119.97,35.88
120.45,36.38
120,36.28
118.05,36.78
117.57,34.86
117.17,35.09
118.49,37.46
118.54,37.59
118.25,37.49
116.29,37.45
116.8,37.64
117.22,37.74
117.15,37.31
117.2,36.97
116.66,36.95
116,36.95
116.58,37.34
117.37,37.37
116.86,37.2
116.76,36.79
116.44,37.16
116.08,37.2
118.03,37.36
117.97,37.47
118.41,37.04
118.12,36.95
117.75,36.89
117.58,37.65
118.14,37.7
118.12,37.12
117.66,37.18
117.51,17.49
117.58,37.73
119.1,36.62
119.22,36.77
119.97,36.77
119.42,35.99
119.2,36.42
118.53,36.5
118.73,36.86
119.41,36.86
119.75,36.38
119.2,35.74
118.83,36.69
118.47,36.69
121.39,37.52
121.59,37.38
122.05,37.2
121.17,36.76
120.71,36.97
120.83,37.28
119.93,37.18
120.73,37.91
122.1,37.5
121.27,37.49
122.41,37.16
121.52,36.89
120.53,36.86
120.38,37.35
120.51,37.64
120.75,37.8
118.35,35.05
118.64,35.78
119.46,35.42
118.73,34.89
118.03,34.84
117.63,35.49
118.17,36.18
118.47,35.54
118.83,35.57
118.83,35.17
118.35,34.61
117.97,35.26
117.95,35.7
117.13,36.18
117.67,36.19
116.76,36.24
116.46,36.29
117.67,35.86
117.76,35.91
116.8,35.76
116.3,35.91
116.59,35.38
116.83,35.54
117.27,35.65
116.65,35
116.34,35.41
116.49,35.71
116.98,35.59
116.97,35.39
117.12,34.8
116.32,35.07
115.43,35.24
115.94,35.59
116.08,35.38
116.07,34.82
115.53,34.83
115.5,35.57
116.1,35.8
115.88,34.97
115.57,35.07
115.08,35.31
115.97,36.45
116.23,36.86
116.23,36.32
115.67,36.24
115.72,36.68
116.27,36.58
115.78,36.11
115.45,35.47
112.53,37.87
112.65,38.05
111.78,38.05
112.33,37.62
113.3,40.12
113.57,37.85
113.08,36.18
114.08,40.42
114.2,39.47
113.1,39.82
112.82,39.52
112.12,39.53
112.33,40.18
113.72,40.38
113.27,39.75
113.68,39.7
113.18,39.58
112.42,39.32
112.67,40.02
112.7,38.38
112.97,39.07
113.32,38.72
111.9,38.37
111.09,38.01
111.17,39.38
112.17,39.1
112.7,38.73
113.28,39.2
112.95,38.5
111.58,38.7
111.82,38.93
111.47,39.45
112.28,39
112.72,37.68
113.37,38.01
113.68,37.62
113.35,37.07
112.53,37.42
112.18,37.2
111.77,36.83
113.17,37.88
113.62,37.79
113.55,37.33
112.97,37.08
112.33,37.36
111.88,37.03
111.13,37.53
111.22,38.47
111.24,37.86
111.62,38.28
112.14,37.55
112.02,37.42
111.75,37.27
111.8,37.12
111.2,36.97
110.83,37
111.17,37.37
110.95,37.95
110.85,37.45
113.02,36.55
113.4,36.56
113.23,35.11
112.88,35.48
112.38,35.84
112.87,36.13
112.32,36.5
113.22,36.33
112.83,36.83
113.43,36.19
113.27,35.78
112.83,35.52
112.15,35.67
112.87,36.32
112.68,36.75
111.5,36.08
111.53,36.63
112.2,36.15
111.9,36.29
111.68,35.73
111.33,35.63
110.65,36.12
110.72,36.47
111.45,35.03
110.64,36.62
111.68,36.25
111.72,36.57
111.83,35.97
111.43,35.86
110.8,35.97
111.07,36.42
110.97,35.03
111.2,35.37
111.63,35.3
110.68,34.71
110.78,35.15
111.22,35.62
110.7,35.58
111.22,35.12
111.58,35.48
111.2,34.12
110.42,34.88
110.83,35.42
110.97,35.6
113.23,23.16
113.19,23.4
114.2,24.09
113.81,23.13
113.55,23.57
114.25,23.75
113.36,22.95
110.35,20.02
116.69,23.39
110.38,21.2
110.88,21.68
113.11,23.05
113.06,22.61
114.07,22.62
113.85,22.58
113.52,22.3
113.62,24.84
113.58,24.68
113.35,25.14
113.73,25.11
114.33,25.14
114.08,24.78
114.13,24.36
113.52,23.86
113.38,24.17
113.01,23.7
112.65,24.48
112.4,24.77
112.07,24.59
112.28,24.77
114.4,23.09
114.4,23.09
114.28,23.18
114.68,23.73
114.48,24.39
114.89,24.45
115.25,24.09
115.18,23.64
114.7,22.97
113.75,23.04
116.1,24.55
116.1,24.55
117.9,24.59
116.18,24.66
116.7,24.34
116.18,23.78
115.75,23.93
115.75,24.15
116.63,23.68
116.8,23.48
116.63,23.68
117.01,23.7
117.03,23.44
116.61,23.27
116.29,23.07
117.64,22.95
117.33,22.98
116.17,23.29
115.82,23.45
116.35,23.55
113.11,23.05
112.89,23.18
113.24,22.84
113.38,22.52
113.25,22.2
113.02,22.52
112.94,22.76
112.68,22.36
112.78,22.27
112.29,22.21
112.76,21.71
110.27,21.63
110.59,21.64
110.83,21.95
110.9,22.36
111.78,22.16
111.95,21.85
110.99,21.52
110.78,21.43
110.17,20.34
110.07,20.91
110.24,21.39
112.44,23.05
112.44,23.05
112.18,23.93
112.43,23.14
112.68,23.36
112.2,22.68
112.02,22.93
111.56,22.77
111.51,23.23
111.75,23.15
111.48,23.45
108.33,22.84
109.4,24.33
110.28,25.29
111.34,23.51
106.75,22.11
108.49,22.74
108.27,23.17
108.2,23.73
108.59,23.44
108.8,23.22
109.2,22.69
107.92,22.65
107.37,22.42
107.08,22.12
106.84,22.36
107.21,22.85
107.12,23.08
107.68,23.18
108.06,24.7
108.26,24.83
108.9,24.79
108.64,24.47
107.36,24.53
107.05,24.55
107.16,25.01
107.54,24.98
108.09,23.94
107.25,24.15
108.89,23.82
109.24,24.67
109.37,24.24
109.74,24.49
109.7,23.98
109.66,23.6
109.34,24.27
109.24,23.76
108.66,24.07
109.24,25.07
109.58,25.8
110.18,24.14
110.22,25.22
110.33,25.42
110.66,25.6
110.66,26.03
111.06,25.96
111.14,25.49
110.81,24.85
110.66,24.64
110.38,24.51
109.98,24.99
110.02,25.78
111.22,23.51
111.3,24.53
110.26,24.83
111.54,24.44
111,22.95
110.9,23.36
110.54,24.22
110.8,24.18
110.14,22.64
110.07,23.38
110.4,23.55
110.53,22.87
110.33,22.71
110.25,22.33
109.98,22.27
109.6,23.11
109.12,21.49
108.61,21.96
109.29,22.44
109.56,22.27
109.2,21.33
107.98,22.16
108.35,21.78
106.62,23.91
106.55,24.35
106.56,24.78
106.9,23.75
107.12,23.62
107.59,23.33
106.6,23.34
106.41,23.15
105.85,23.42
105.08,24.51
106.24,24.31
105.34,24.8
87.68,43.77
84.77,45.59
85.94,44.27
89.19,42.91
88.63,42.77
90.25,42.82
93.44,42.78
94.65,43.28
93,43.6
86.06,41.68
86.35,42.31
86.84,42.23
86.53,41.95
86.24,41.36
84.25,41.77
86.55,42.05
79.94,37.12
82.63,37.07
80.78,37.04
81.63,36.86
80.17,37.12
78.29,37.06
79.74,37.31
80.29,41.15
80.24,41.29
81.84,41.82
82.97,41.68
82.63,41.55
82.9,41.25
80.34,40.64
79.06,40.55
79.25,41.22
75.94,39.52
78.59,39.78
76.78,39.46
76.67,39.23
77.62,38.95
77.25,38.45
77.26,38.2
77.42,37.89
76.05,39.41
76.17,38.91
75.83,39.42
75.22,37.76
76.12,39.73
78.42,41.91
75.94,39.14
75.18,39.7
87.31,44.05
87.94,44.14
89.52,44.02
89.14,44
87.68,43.97
86.22,44.28
86.92,44.18
90.34,43.8
82.1,44.93
82.92,44.67
81.08,44.95
81.33,43.91
82.53,43.82
83.27,43.41
82.23,43.35
84.89,44.45
81.81,43.23
81.08,43.15
80.87,44.07
81.12,43.82
82.96,46.74
83.62,46.52
84.62,44.45
83.59,45.92
82.94,46.21
85.56,44.29
85.13,46.78
88.14,47.86
90.37,46.71
89.44,47.05
87.51,47.15
85.84,47.42
86.92,47.7
86.41,48.05
118.78,32.04
118.83,31.95
118.83,32.36
118.62,32.07
117.2,34.26
119.16,34.59
120.86,32.01
120.62,31.32
120.29,31.59
119.95,31.79
116.57,34.79
116.93,34.73
119.11,34.83
118.75,34.54
118.33,34.38
117.97,34.3
117.94,33.89
117.2,34.26
119.02,33.59
119.23,34.3
119.36,34.09
118.79,34.12
118.3,33.96
118.68,33.73
118.05,33
119.26,33.77
119.02,33.62
119.15,33.5
118.85,33.28
118.23,33.46
119.02,33.01
120.13,33.38
119.84,34.01
119.79,33.78
120.26,33.77
119.77,33.46
119.56,34.2
120.45,33.19
120.31,32.84
120.45,32.57
120.56,32.39
121.18,32.33
121.66,31.8
121.15,31.89
120.86,32.01
119.42,32.39
119.32,33.23
119.82,32.93
119.45,32.78
120.02,32.16
120.15,32.51
119.9,32.49
120.26,32.03
119.55,32.43
119.42,32.39
119.16,32.27
119.44,32.2
119.44,32.2
119.81,32.24
119.55,32
119.95,31.78
119.82,31.36
119.56,31.74
119.48,31.43
119.16,31.95
119.02,31.65
118.87,31.32
120.26,31.91
120.55,31.86
120.74,31.64
121.1,31.45
120.95,31.39
120.62,31.32
120.63,31.16
115.89,28.68
115.8,28.69
117.22,29.3
113.85,27.6
115.97,29.71
116.56,29.9
116.23,29.75
116.19,29.29
116.03,29.47
115.82,29.04
115.75,29.33
115.65,29.68
115.09,29.26
114.55,29.04
117.97,28.47
117.83,29.25
117.58,28.96
118.25,28.68
118.2,28.45
117.71,28.32
117.62,28.42
117.02,28.23
117.2,28.3
116.82,28.22
117.08,28.7
117.12,28.97
116.68,29
116.69,28.7
117.43,28.42
114.38,27.81
114.44,28.11
114.37,28.53
114.78,28.4
114.91,28.25
115.55,28.86
115.38,28.71
115.38,28.42
115.7,28.19
115.54,28.07
114.92,27.81
114.68,27.82
115.37,28.88
116.34,28
116.29,27.95
116.77,27.92
117.06,27.7
116.91,27.3
116.52,27.22
116.62,27.56
116.2,27.55
116.05,27.75
115.82,27.44
116.61,28.23
116.26,28.37
114.97,27.12
115.4,27.77
115.15,27.56
115.14,27.22
115.42,27.33
114.88,26.81
114.77,26.47
114.5,26.33
113.97,26.71
114.23,26.96
113.94,27.14
114.62,27.39
114.17,26.57
114.92,25.85
116.32,26.84
116.32,26.34
116,26.46
115.33,26.32
115.39,25.96
116.02,25.89
115.79,25.58
115.41,25.15
115.64,24.96
115.02,24.7
114.79,24.91
114.53,24.76
114.94,25.39
114.02,25.85
114.75,25.66
114.55,25.8
114.31,25.69
114.36,25.39
114.48,38.03
118.02,39.63
114.54,38.42
114.38,38.31
115.18,37.94
115.03,38.03
114.84,38.03
114.58,37.62
114.78,37.76
114.13,38.03
114.03,38.08
114.67,38.33
114.56,38.13
115.2,38.2
114.96,38.16
114.35,37.65
114.5,37.74
114.64,38.87
114.24,38.2
114.47,36.6
114.5,36.77
114.92,36.78
115.4,36.47
114.94,36.37
114.68,36.43
115.14,36.28
113.67,36.57
113.85,36.95
115.18,36.84
114.94,36.49
114.8,36.56
114.62,36.35
114.37,36.37
114.2,36.7
114.48,37.05
114.68,37.49
114.9,37.62
114.75,37.35
115.5,36.87
115.37,37.37
115.03,37.22
114.68,37.11
114.52,36.94
114.5,37.43
114.5,37.28
115.22,37.53
115.67,37.07
115.08,36.97
115.14,37.06
115.02,37.06
114.71,37
115.48,38.85
115.71,39.39
115.98,39.48
115.78,39.28
115.86,39.06
115.92,38.92
115.58,38.49
115.46,38.46
114.02,38.52
114.18,38.85
114.97,38.75
114.67,39.37
115.49,39.35
115.84,39.34
116.1,38.98
115.65,39.02
115.78,38.68
115.3,38.41
115.47,38.76
115.14,38.71
114.68,38.62
115.12,38.84
115.45,38.95
114.87,40.82
114.6,41.87
115.82,40.92
115.54,40.4
114.53,39.83
115.03,40.63
114.7,41.15
115.68,41.68
115.25,40.98
115.2,40.37
114.15,40.12
114.38,40.67
113.95,41.05
114.73,40.84
117.93,40.97
117.72,41.95
118.68,41.02
118.47,40.62
117.48,40.42
117.53,40.95
117.7,41.32
118.93,40.43
116.63,41.2
119.57,39.95
118.3,40.15
118.69,40.02
119.15,39.72
118.85,39.89
118.67,39.49
117.9,39.9
118.54,39.31
117.97,40.2
119.22,39.88
118.9,39.43
118.73,39.74
118.1,39.58
118.13,39.82
116.7,39.53
116.69,39.52
117.06,39.97
117,39.76
116.38,39.12
116.29,39.44
116.63,38.7
116.45,38.87
116.48,39.32
116.98,39.98
116.83,38.33
117.33,38.4
117.22,38.07
116.37,37.65
116.52,37.89
115.82,38.43
116.07,38.45
116.56,38.08
116.27,38.02
116.8,38.58
117.85,38.17
116.7,38.05
116.08,38.72
116.12,38.2
117.1,38.06
115.72,37.72
115.74,38.24
116.14,37.87
116.26,37.69
115.72,37.52
115.56,38.02
115.5,38.22
115.96,38.03
115.9,37.81
115.96,37.36
115.56,37.59
113.65,34.76
113.35,34.79
114.35,34.79
113.29,33.75
112.44,34.7
113.21,35.24
114.17,35.9
114.77,34.56
114.17,34.41
113.71,34.4
113.02,34.46
114.46,34.48
114,34.73
113.35,34.51
112.96,34.76
114.81,34.69
113.85,35.31
114.05,35.44
114.04,35.03
113.63,35.27
113.06,34.94
112.57,35.08
113.05,35.16
113.77,35.46
114.19,35.14
113.96,35.05
113.38,35.1
112.77,34.92
112.92,35.08
113.42,35.24
114.35,36.1
115.21,36.08
115.46,35.9
115.83,36
114.52,35.57
114.54,35.67
114.17,35.6
114.88,35.95
115.1,35.89
114.98,35.71
114.67,35.19
114.35,35.92
113.81,36.06
115.65,34.44
116.13,34.22
115.29,34.08
115.04,34.46
115.87,34.4
116.37,33.94
115.31,34.44
115.13,34.65
114.63,33.63
114.59,33.54
114.38,34.05
115.48,33.86
114.88,33.74
115.06,33.41
114.5,33.79
114.85,34.06
115.17,33.63
114.9,33.44
113.81,34.02
114.17,34.11
113.98,33.6
113.46,33.86
112.88,33.74
113.19,33.98
114.02,33.56
113.77,34.22
113.94,33.81
113.58,33.44
113.35,33.62
113.04,33.86
113.47,34.16
114.02,32.98
114.02,32.83
114,33.38
114.35,33
114.97,32.75
113.31,32.72
113.98,33.15
114.26,33.25
114.62,32.97
114.38,32.62
114.08,32.13
114.72,32.35
115.68,32.17
115.04,32.13
114.83,31.62
114.53,32.21
115.41,32.44
115.42,31.81
114.91,32.02
112.53,33.01
112.98,33.25
112.83,32.7
112.36,32.51
112.08,32.68
111.47,33.14
112.4,33.49
112.92,33.05
113.4,32.37
112.23,33.03
111.83,33.05
111.5,33.31
111.19,34.76
112.42,34.84
112.83,34.17
112.46,34.16
112.07,34.14
111.6,33.81
110.85,34.52
111.75,34.76
111.92,34.73
112.77,34.73
112.42,34.43
112.15,34.51
111.65,34.39
111.03,34.06
111.19,34.76
112.14,34.75
120.19,30.26
120.3,30.43
119.95,30.07
119.27,29.49
119.72,30.23
120.25,30.16
119.64,29.8
119.05,29.61
121.56,29.86
121.72,29.96
120.65,28.01
120.65,28.01
120.68,28.16
121.12,27.84
120.55,27.68
119.7,27.57
120.94,28.14
120.62,27.8
120.08,27.08
120.36,27.53
120.1,30.86
121.02,30.7
120.54,30.64
119.68,30.68
120.92,30.84
120.76,30.77
120.92,30.53
120.69,30.53
120.08,30.54
119.91,30.01
122.11,30.03
122.2,30.26
122.45,30.72
122.3,29.97
121.56,29.86
121.8,29.48
121.41,29.66
121.23,30.18
121.42,29.3
121.16,30.04
120.58,30.01
120.89,29.49
120.23,29.71
120.87,30.03
120.81,29.6
121.44,28.67
121.13,28.8
121.38,29.11
121.36,28.36
120.73,28.85
121.03,29.15
121.27,28.64
121.23,28.14
119.92,28.45
120.28,28.45
119.06,27.61
119.25,28.59
120.6,28.66
119.56,28.12
119.13,28.08
119.48,28.46
119.64,29.12
119.88,29.46
120.23,29.27
119.81,28.9
118.61,28.74
118.39,29.15
118.88,28.97
119.48,29.19
120.06,29.32
120.02,28.92
118.5,28.9
110.35,20.02
110.33,19.98
110.72,19.61
110.31,19.68
110.46,19.25
110.39,18.8
110.1,19.36
110,19.75
109.57,19.52
109.69,19.91
109.7,18.64
109.44,19.23
109.83,19.05
110.02,18.48
109.5,18.25
109.17,18.73
108.64,19.09
109.03,19.25
114.1,22.2
113.33,22.13
121.5,25.05
120.37,22.64
121.73,25.14
120.67,24.15
120.19,22.98
121.75,24.75
121.3,25
120.96,24.81
114.31,30.52
114.33,30.35
114.02,30.57
115.09,30.2
110.79,32.65
112.24,30.32
111.3,30.7
112.14,30.02
113.91,31.92
114.36,30.88
113.59,30.63
113.73,31.02
113.81,31.62
114.09,31.56
113.6,30.94
113.69,31.25
114.87,30.38
114.87,30.44
114.8,31.84
114.61,31.29
115,31.17
115.37,30.79
115.22,30.46
115.3,30.24
115.93,30.09
115.56,29.85
115.57,30.75
114.28,29.87
115.22,29.83
114.52,29.6
113.8,29.23
113.91,29.97
114.04,29.54
113.85,29.71
112.19,31.02
112.18,30.35
112.58,31.17
113.11,31.03
112.9,29.83
112.41,29.73
113,28.21
112.8,28.37
113.16,27.83
112.91,27.87
112.61,26.89
111.5,27.22
113.09,29.37
113.42,29.48
113.56,29.71
113.05,28.8
112.87,28.68
112.55,29.52
113.63,28.16
113.5,27.67
113.32,27.01
113.54,26.79
113.77,26.49
112.5,27.75
113,25.79
113,25.79
113.27,26.71
113.11,26.13
113.39,25.95
113.91,25.08
113.68,25.54
113.96,25.41
112.55,25.27
112.35,25.56
112.72,25.73
112.84,26.41
112.61,26.89
112.86,27.25
112.95,27.1
112.39,26.38
111.85,26.59
112.14,26.8
112.61,26.89
111.63,26.22
111.63,26.22
112.21,25.91
111.95,25.6
112.16,25.37
111.64,25.96
111.33,25.41
111.57,25.52
111.28,26.41
111.79,24.97
110.84,26.44
110.61,26.73
111.04,27.13
110.14,25.59
110.57,27.06
110.3,26.37
111.96,27.71
111.66,27.68
111.46,27.33
112.18,27.44
111.41,27.68
111.73,27.25
111.29,27.73
109.95,27.52
110.14,27.33
110.18,28.02
110.39,28.46
110.57,27.92
109.71,26.86
109.68,26.57
109.96,27.1
109.78,27.44
109.79,27.87
109.77,26.16
109.16,27.37
109.71,28.3
109.84,29
110.16,29.38
110.48,29.13
109.91,28.62
110.73,28.29
109.43,27.92
109.46,28.59
109.64,28.7
109.42,29.64
111.69,29.05
111.64,29.44
111.75,29.65
112.16,29.41
111.87,29.64
111.97,28.9
111.47,28.9
111.09,29.41
111.35,29.59
112.33,28.6
112.39,29.37
112.36,28.83
112.55,28.27
111.2,28.38
112.11,28.51
103.73,36.03
103.25,36.73
104.09,35.87
101.94,38.23
103.97,36.32
104.57,35.57
105.08,35.72
104.61,34.98
103.88,35.39
104.71,36.54
105.27,35.24
104.19,35.17
106.68,35.51
107.61,35.1
106.65,35.21
105.73,35.51
107.38,35.31
107.05,35.27
106.06,35.2
107.88,36.03
108,36.44
108.43,35.5
107.22,35.7
107.33,36.57
108.02,35.81
107.94,35.17
105.69,34.6
106.11,33.78
105.15,34.22
104.88,34.69
105.69,34.89
106.12,34.73
106.28,33.9
105.28,34.02
105.35,34.7
104.48,34.87
106.23,35
104.94,33.43
104.38,34.06
105.58,33.33
105.7,33.75
104.7,32.95
103.35,34.69
104.38,33.81
102.04,33.97
102.46,35.21
103.54,34.61
103.23,34.08
102.5,34.6
103.22,35.62
103.34,35.97
103.31,35.43
103.68,35.39
103.54,35.46
103.39,35.68
104.04,34.41
102.85,35.74
102.61,37.94
103.08,38.62
102.86,37.43
104.05,37.14
102.84,37.24
100.46,38.93
100.85,38.43
100.17,39.14
101.19,38.79
99.84,39.14
99.57,38.86
97.58,39.81
98.5,39.71
94.71,40.13
98.92,39.97
95.77,40.51
94.25,38.46
94.89,39.49
119.3,26.08
119.14,26.16
118.1,24.46
118.15,24.74
118.16,26.65
118.16,26.65
118.32,27.05
118.55,27.92
117.48,27.34
117.8,26.8
118.02,27.76
117.34,27.54
118.77,27.53
118.85,27.38
119.52,26.65
119.65,27.09
119.53,26.2
120.2,27.34
120,26.89
118.74,26.59
119.55,26.49
119.5,27.47
119.36,27.12
118.98,26.92
119.89,27.25
119,25.44
118.7,25.37
119.39,25.73
119.52,25.96
118.95,25.88
119.78,25.51
118.86,26.21
118.58,24.93
118.57,24.82
118.39,24.96
118.78,25.04
118.18,25.07
118.3,25.34
118.24,25.5
118.34,24.43
117.35,24.52
117.79,24.44
117.61,24.12
117.16,23.73
117.3,24.38
117.34,23.99
117.35,24.51
117.75,24.62
117.4,23.72
117.53,25
117.01,25.12
116.41,25.43
116.81,24.76
116.37,25.85
116.1,25.11
116.75,25.72
117.4,25.3
117.61,26.23
118.17,26.18
116.64,26.26
117.83,25.69
117.37,25.97
117.77,26.41
117.45,26.73
116.81,26.12
116.82,26.85
117.15,26.92
117.18,26.36
91.11,29.97
91.24,30.2
91.05,30.51
91.77,29.77
90.14,29.44
94.13,29.18
95.26,29.22
91.39,29.63
90.7,29.39
90.96,29.67
94.25,29.59
93.25,29.92
92.1,31.47
94.1,31.96
93.68,31.53
90.05,31.35
93.46,30.63
92.3,31.08
93.71,31.92
91.68,32.29
88.7,30.94
97.14,31.18
98.29,30.86
97.9,29.68
97.49,28.62
95.76,30.81
95.63,31.42
95.75,29.92
89.19,31.53
97.56,30.69
98.68,29.64
96.95,30.04
94.69,30.94
96.57,31.2
91.76,29.18
92.6,29.09
92.11,29.08
91.91,27.98
91.65,29.04
90.96,29.25
90.33,29.96
92,29.26
93.11,29.06
92.42,28.46
91.4,28.49
90.83,28.42
91.26,29.22
88.82,29.28
87.77,28.38
87.62,29.1
85.94,28.19
88.25,29.43
84.15,29.66
89.67,28.57
88.93,27.55
88.5,28.29
89.02,29.71
88,28.87
87.11,28.57
85.29,28.94
87.22,29.3
89.63,28.94
89.77,29.21
89.16,29.11
85.3,29.38
80,32.08
81.13,32.45
79.76,31.47
85.16,31.06
79.61,33.44
84.1,32.33
81.18,30.37
106.71,26.57
104.82,26.58
104.82,26.58
104.64,25.81
105.47,26.21
106.9,27.7
107.19,27.95
107.6,28.89
107.72,27.97
107.88,27.22
105.69,28.57
106.8,28.16
107.43,28.56
107.87,28.54
107.5,27.76
106.41,27.81
106.2,28.33
109.21,27.73
108.91,27.24
108.23,27.94
108.13,28.27
109.2,27.52
108.82,27.68
108.24,27.52
108.41,28.02
108.48,28.57
109.18,28.17
105.29,27.32
106.04,27.03
105.76,26.66
104.71,27.13
105.61,27.16
106.22,27.46
105.38,26.77
104.28,26.87
105.92,26.25
106.73,27.1
106.46,26.56
105.75,26.32
106.95,27.06
106.59,26.84
106.26,26.42
105.75,26.08
106.06,25.75
105.62,25.94
104.91,25.1
104.96,25.79
105.63,25.39
106.09,25.17
105.79,25
105.49,25.11
105.18,25.44
105.21,25.83
107.97,26.59
108.11,27.03
108.41,27.06
109.2,26.89
108.58,26.64
109.14,26.24
108.9,25.76
107.58,26.49
107.89,26.89
108.68,26.98
108.72,27.21
109.18,26.7
108.32,26.68
108.5,25.94
108.07,26.38
107.79,26.21
107.53,26.72
107.22,26.58
107.48,27.08
107.55,25.83
106.45,26.03
106.66,26.14
107.88,25.42
107.51,26.7
107.54,25.84
106.74,25.43
106.98,26.46
107.86,26
123.38,41.8
122.83,42
122.7,41.52
121.62,38.92
121.7,39.13
121.97,39.63
121.95,39.55
122.58,39.28
122.85,41.12
122.75,40.85
122.4,41.4
123.97,41.97
125.02,41.72
124.9,42.13
123.73,41.3
125.33,41.28
121.15,41.13
121.35,41.17
121.22,41.55
122.12,41.7
121.8,41.6
120.83,40.77
120.68,40.63
120.32,40.35
124.37,40.13
124.13,39.97
123.25,40.3
124.05,40.47
124.77,40.75
121.65,42
122.52,42.42
122.18,40.65
122.37,40.42
122.03,41.02
122.06,41
123.17,41.28
123.34,41.43
123.85,42.32
124.03,42.53
124.13,42.8
123.5,42.48
123.33,42.75
123.37,42.52
124.7,42.77
120.42,41.58
119.78,40.82
120.75,41.82
119.37,41.27
119.63,41.38
106.54,29.59
106.56,29.41
106.64,29.01
107.04,29.86
106.28,29.26
106.22,30.03
105.8,30.16
106.03,29.86
106.21,29.62
105.59,29.4
105.71,29.75
105.91,29.38
108.95,34.27
108.97,34.18
109.11,35.09
108.98,34.91
107.15,34.38
107.39,34.53
107.13,34.65
106.86,34.91
107.8,34.69
107.63,34.46
107.87,34.38
108.22,34.28
107.76,34.29
107.3,34.09
106.51,33.93
109.77,38.3
110.51,38.83
111.07,39.05
110.48,38.04
110.23,37.78
110.73,37.49
110.24,37.49
110.15,37.11
110.05,37.45
109.32,37.97
108.79,37.61
107.59,37.6
109.47,36.6
109.34,36.88
109.65,37.16
110.18,36.87
110.02,36.59
110.15,36.04
109.86,35.6
109.42,35.76
109.11,35.43
109.27,35.6
109.37,36
109.37,36.29
108.78,36.84
108.22,36.93
108.72,34.36
108.43,34.5
108.14,34.71
108.09,35.04
107.8,35.22
108.33,35.13
108.57,34.81
108.84,34.53
108.94,34.62
109.1,34.55
108.61,34.12
108.22,34.18
108.49,34.32
108.25,34.54
109.5,34.52
109.59,34.97
109.6,35.18
109.93,35.2
110.45,35.47
110.15,35.24
109.96,34.82
110.25,34.56
110.09,34.58
109.77,34.53
109.32,34.17
109.22,34.38
109.17,34.76
109.96,33.88
110.15,34.11
110.35,33.71
110.88,33.54
109.91,33.55
109.16,33.45
109.14,33.69
109.02,32.7
109.35,32.83
110.06,32.83
109.37,32.41
109.51,31.91
108.89,32.3
108.55,32.56
108.53,32.9
108.26,33.05
108.33,33.34
108.04,33.07
106.95,33.65
107.32,33.16
107.56,33.23
108,33.55
107.77,33
107.91,32.56
106.93,33
106.25,32.82
106.68,33.16
106.16,33.34
101.74,36.56
101.67,36.92
102.09,36.47
101.57,36.49
102.38,36.49
102.8,36.3
101.28,36.72
101.95,36.84
102.3,36.11
102.46,35.84
101.62,37.37
100.99,36.89
100.17,37.32
100.22,38.2
102,35.54
102,35.92
101.5,35.03
101.62,34.75
100.61,36.27
101.47,36.02
100.75,35.57
100.63,35.24
99.99,35.6
100.26,34.49
99.89,33.95
101.47,33.46
100.73,32.92
99.68,33.74
98.26,34.92
96.97,33.03
97.12,33.35
96.47,32.23
95.3,32.92
95.6,33.86
95.5,34.52
94.9,36.41
98.46,36.9
98.13,36.3
99.03,37.28
126.63,45.75
123.97,47.33
130.3,47.33
131.17,46.65
130.97,45.3
125.03,46.58
128.92,47.73
130,48.93
128.08,47.98
127,46.63
127.12,47.22
126.97,47.47
127.5,46.87
126.3,46.28
125.98,46.07
125.25,45.72
125.07,45.53
125.33,46.42
125.88,47.18
126.13,46.68
126.5,46.83
127.53,50.22
127.53,50.22
126.17,48.5
126.8,49.76
126.5,48.22
127.5,49.22
128.42,49.57
125.2,49.17
130.35,46.83
130.68,47.02
130.83,47.58
131.83,47.3
132.02,47.23
132.5,47.67
134.15,48.33
134,46.78
130.83,45.82
132.17,46.33
131.13,46.7
130.53,45.75
130.53,46.25
129.55,46.33
129.92,46.73
129.58,44.6
130.23,45.3
131.04,45.27
131.85,45.53
133.97,45.75
131.17,44.38
131.12,44.07
130.5,44.9
129.47,44.35
129.35,44.57
126.95,45.52
126.58,46
127.38,46.08
127.48,45.75
128.03,45.95
128.7,45.98
128.8,45.83
128.35,45.47
127.95,45.22
127.17,44.93
126.32,45.53
124.4,47.8
124.85,48.48
125.87,48.03
126.22,48.03
126.07,47.62
125.3,47.92
124.87,47.18
123.45,46.4
123.18,47.35
123.48,47.9
124.44,46.86
124.07,50.42
126.6,51.72
124.7,52.32
122.37,53.48
[/mdx_fold]

分类
程序和算法

大数据与机器学习 基础篇 回归

回归,Regression,是一种归纳的思想——当看到大量的事实所呈现的样态,推断出原因是如的,当看到大量的数字对(pair)是某种样态,推断出它们之间蕴含的关系是如何。

注:本文中用到的Python及其模块安装教程参见


 线性回归

线性回归是利用数理统计学中的回归分析来确定两种或两种以上变量间相互依赖的定量关系的一种统计分析方法。其表达形式如下:

\(y=ax+b+e\)

(e为误差服从均值为0的正态分布)

假设有多个x和y样本值,尝试用y=ax+b+e来拟合,可以得到:

\(\vert e\vert=\vert ax+b-y\vert\)

误差大小其实是猜想的ax+b的值和观测到的y的值之间的差值。试着把所有的|e|都求和:

\(Q=\sum_{i=1}^n(ax_i+b-y_i)^2\)

现在的问题是让a和b分别等于什么值时Q最小,即:

\(\frac{∂Q}{∂a}=0\)且\(\frac{∂Q}{∂b}=0\)

分别求导后:

\(
\left\{
\begin{array}{c}
\frac{∂Q}{∂a}=2\sum_{i=1}^n[x_i(ax_i+b-y_i)]=0\\
\\
\\
\frac{∂Q}{∂b}=2\sum_{i=1}^n(ax_i+b-y_i)=0
\end{array}
\right.
\)

推导出:

\(
\left\{
\begin{array}{c}
\sum_{i=1}^nx_i^2+b\sum_{i=1}^nx_i=\sum_{i=1}^ny_ix_i(1)\\
\\
\\
a\sum_{i=1}^nx_i+nb=\sum_{i=1}^ny_i(2)
\end{array}
\right.
\)

把(2)式做如下处理:

\((2)·\frac{\sum_{i=1}^nx_i}{n}-(1)\)

得到下面的等式:

\(a\Biggl( \sum_{i=1}^nx_i·\frac{\sum_{i=1}^nx_i}{n}-\sum_{i=1}^nx_i^2\Biggr)=\sum_{i=1}^ny_i·\frac{\sum_{i=1}^nx_i}{n}-\sum_{i=1}^ny_ix_i\)

推导出来a和b的表达式:

\(
\left\{
\begin{array}{c}
a=\frac{\frac{\sum_{i=1}^ny_i·\sum_{i=1}^nx_i}{n}-\sum_{i=1}^ny_ix_i}{\frac{\sum_{i=1}^nx_i·\sum_{i=1}^nx_i}{n}-\sum_{i=1}^nx_i^2}\\
\\
\\
b=\frac{\sum_{i=1}^ny_i-a\sum_{i=1}^nx_i}{n}
\end{array}
\right.
\)

其中,a就是斜率,b是截距。

在Python中,可以直接使用最小二乘法进行线性拟合。例如:

[mdx_table]
x
1
2
3
4
5
6
7
8
9
—–
y
0.199
0.389
0.580
0.783
0.980
1.177
1.380
1.575
1.771
[/mdx_table]

import numpy as np
import matplotlib.pyplot as plt

x=[1,2,3,4,5,6,7,8,9]
y=[0.199,0.389,0.580,0.783,0.980,1.177,1.380,1.575,1.771]

t1=t2=t3=t4=0
n=len(x)
for i in range(n):
    t1+=y[i]#\sum y
    t2+=x[i]#\sum x
    t3+=x[i]*y[i]#\sum xy
    t4+=x[i]**2#\sum x^2

a=(t1*t2/n-t3)/(t2*t2/n-t4)
b=(t1-a*t2)/n

x=np.array(x)
y=np.array(y)

plt.plot(x,y,'o',label='Original data',markersize=10)
plt.plot(x,a*x+b,'r',label='Fitted line')
plt.show()
  • 拟合的效果如图:

也可以直接用numpy包中的内置函数,效果与上图相同。

import numpy as np
import matplotlib.pyplot as plt

x=[1,2,3,4,5,6,7,8,9]
y=[0.199,0.389,0.580,0.783,0.980,1.177,1.380,1.575,1.771]

A=np.vstack([x,np.ones(len(x))]).T

#调用最小二乘法函数
a,b=np.linalg.lstsq(A,y)[0]

x=np.array(x)
y=np.array(y)

plt.plot(x,y,'o',label='Original data',markersize=10)
plt.plot(x,a*x+b,'r',label='Fitted line')
plt.show()

非线性回归

有些非线性回归模型,经过适当的数学变换,便能得到它的线性化表达形式,但对另外一些非线性回归模型,仅仅做变量变换根本无济于事。属于前一种情况的非线性回归模型一般称为内蕴的线性回归,而后者则称为内蕴的非线性回归。

这里只讨论前者这种简单的情况,例如:

  • 表 20世纪60年代世界人口状况

[mdx_table header=”true”]
年份
世界人口(亿)
—–
1960
29.72
—–
1961
30.61
—–
1962
31.51
—–
1963
32.13
—–
1964
32.34
—–
1965
32.85
—–
1966
33.56
—–
1967
34.20
—–
1968
34.83
[/mdx_table]

根据马尔萨斯人口模型(其中s是人口,t是年份,e是自然常数):

\(s=α·e^{βt}\)

对等式两边同时取ln

\(lns=βt+lnα\)

这里实际上用的还是线性回归模型,相当于:

\(y=ax+b\)

\(
\left\{
\begin{array}{c}
lns=y\\
\\
β=a\\
\\
lnα=b
\end{array}
\right.
\)

直接调用最小二乘法拟合:

import numpy as np
import matplotlib.pyplot as plt

x=[1960,1961,1962,1963,1964,1965,1966,1967,1968]
y=[29.72,30.61,31.51,32.13,32.34,32.85,33.56,34.20,34.83]

logy=np.log(y)

A=np.vstack([x,np.ones(len(x))]).T

#调用最小二乘法函数
a,b=np.linalg.lstsq(A,logy)[0]

x=np.array(x)

plt.plot(x,logy,'o',label='Original data',markersize=10)
plt.plot(x,a*x+b,'r',label='Fitted line')
plt.show()
  • 线性拟合的效果如图:

求出β和α

beta=a
alpha=pow(np.e,b)

y=np.array(y)

plt.plot(x,y,'o',label='Original data',markersize=10)
plt.plot(x,alpha*pow(np.e,beta*x),'r',label='Fitted line')
plt.show()

plt.axis([1960,2010,25,100])
plt.plot(x,y,'o',label='Original data',markersize=10)
x=np.arange(1960,2010)
plt.plot(x,alpha*pow(np.e,beta*x),'r',label='Fitted line')
plt.show()
  • 非线性拟合的效果如图:

  • 预测的效果如图:


想了解更多关于大数据和机器学习:

分类
程序和算法

大数据与机器学习 入门篇

注:本文中用到的Python及其模块安装教程参见


大数据产业概述

数据生命周期中的环节

什么是数据

数据是承载一定的信息的符号。

什么是信息?[1]

信息是用来消除随机不定性的东西。


数学基础:统计与分布

加和值

\(\sum_{i=0}^n Xi\)

平均值

\(\overline X=\frac{ \sum_{i=0}^n Xi } {n}\)

标准差

\(σ=\sqrt{\frac{1}{n}\sum_{i=0}^n (Xi-\overline X)^2}\)

加权平均

\(\overline X=\frac{ \sum_{i=0}^n Xi*f(Xi) } {\sum_{i=0}^n f(Xi)}\)

欧式距离

\(d=\sqrt{\sum_{i=0}^n (Xi1-Xi2)^2}\)

曼哈顿距离

\(d=\sum_{i=0}^n \vert Xi1-Xi2\vert\)

同比和环比

同比:相邻大周期的相同小周期的比较。

环比:相邻小周期的比较。

抽样

抽样(Sampling)是一种非常好的了解大量样本空间分布情况的方法,样本越大则抽样带来的成本减少的收益就越明显。

抽样对象要更加具有代表性和分散性,这样才会体现出与整个样本空间更为相近的分布特点。

高斯分布

概率函数:\(f(x)=\frac{1}{\sqrt{2π}σ}exp(-\frac{(x-μ)^2}{2σ^2})\)

  • X的分布:
  • (μ-σ , μ+σ): 68.2%
  • (μ-2σ , μ+2σ): 95.4%
  • (μ-3σ , μ+3σ): 99.6%

泊松分布

概率密度函数:\(P(X=k)=\frac{λ^k}{k!}e^{-λ}\)

参数λ是单位时间(或单位面积)内随机事件的平均发生率。泊松分布适合于描述单位时间内随机事件发生的次数。

  • 泊松分布适用的事件需要满足以下3个条件:
    1. 这个事件是一个小概率事件。
    2. 事件的每次发生是独立的不会相互影响。
    3. 事件的概率是稳定的。

例子:

已知有一个书店,售卖许多图书,其中工具书销售一直较稳定且数量较少(概率较小的事件),新华字典平均每周卖出4套。作为书店老板,新华字典应该备多少本为宜?

每周卖出的新华字典数量k满足λ为4的泊松分布:

  • 表1:不同k值对应的累计概率

[mdx_table header=”true”]
k值
概率
累积概率
—–
1
7.33%
733%
—–
2
14.7%
22.03%
—–
3
19.5%
41.53%
—–
4
19.5%
61.03%
—–
5
15.6%
76.63%
—–
6
10.4%
87.03%
—–
7
5.95%
92.98%
—–
8
2.98%
95.96%
—–
9
1.32%
92.78%
[/mdx_table]

  • 图1:不同k值对应的概率散点图

在泊松分布的例子里,可以看到一个现象:就是k每增加1,在k小于λ的时候,累积函数的增加是很快的,而且每次增加的量比上一次增加的要多;而在k越过λ之后,虽然还在增加,但是每次增加的量比上一次增加的要少,然后越来越少。

伯努利分布

概率函数:\(P(X=k)=C_n^k ·p^k(1-p)^{n-k}\)


指标

指标就是制定的标准,就是为了描述一些对象的状态而制定出来的标准。

  • 指标的选择:
    1. 数字化
    2. 易衡量
    3. 意义清晰
    4. 周期适当
    5. 尽量客观

信息论

信息的定义

首先引用最被大家广泛认可的一种科学性的信息定义——“信息是被消除的不确定性。”[2]

例子:

抛一枚硬币。假设不会出现硬币立在地面上的情况。
结果A说:“硬币落地后正面朝上。”
然后B说:“硬币落地后朝上的面不是反面。”

在我们不知道硬币落地的结果之前,正面朝上的反面朝上的可能性都是存在的,当A告诉我准确的信息之后,那么硬币反面朝上的结果就不存在了,这里“硬币落地后正面朝上”就是信息;而当随机不确定性被消除之后,再被告知的这些信息里就没有消除随机不确定性的因素了,如B说的“硬币落地后朝上的面不是反面”就不是信息。

但如果C说:“这枚硬币最后落在了桌子上”,那么它又是信息,因为它消除了其他的不确定性。

信息量

在信息论中,对信息量是有确定解释并且可以量化计算的,这里的信息量就是一种信息数量化度量的规则。

一段文字有多少信息的想法最早还是在1928年由哈特莱(R.V.L.Hartley)首先提出,他将信息数的对数定义为信息量。

若信源有m种信息,且每个信息是以相等可能产生的,则该信源的信息量可表示如下:

\(I=log_2m\)

如上面提到的抛硬币的例子,因为硬币落地有正面和反面两种可能性,所以m=2,信息量\(I=log_22=1\)。极端情况是,只有一个可能值的时候信息量为0,也就是无须告知也知道结果,即便告知了结果,信息量也为0,如一般情况下硬币抛出后必然会落地,所以“硬币落地”这句话的信息量就是0。

在概率不等的情况下,事件出现的概率越小,信息量越大。Xi表示一个发生的事件,Pi表示这个事件发生的先验概率,则这个事件的信息量为:

\(H(X_i)=-log_2P_i\)

还是上面提到的抛硬币的例子,假设硬币被动过手脚,正面朝上的概率为\(\frac{1}{8}\),反面朝上的概率为\(\frac{7}{8}\),则抛一次硬币之后,正面朝上的信息量为:

\(H(X_i)=-log_2\frac{1}{8}=3\)

反面朝上的信息量为:

\(H(X_i)=-log_2\frac{7}{8}=0.193\)

信息熵

信息熵是信息的杂乱程度的量化描述,公式如下:

\(H(x)=-\sum_{i=1}^np(x_i)log_2P(x_i)\)

信息越确定,越单一,信息熵越小;信息越不确定,越混乱,信息熵越大。

如上面抛硬币的例子中,

  1. 硬币还没有被动过手脚,两面朝上的概率都是\(\frac{1}{2}\):
    信息熵为\(\frac{1}{2}·-log_2\frac{1}{2}+\frac{1}{2}·-log_2\frac{1}{2}=1\)
  2. 硬币已经被动过手脚,正面朝上的概率为\(\frac{1}{8}\),反面朝上的概率为\(\frac{7}{8}\):
    信息熵为\(\frac{1}{8}·-log_2\frac{1}{8}+\frac{7}{8}·-log_2\frac{7}{8}=0.544\)

即知道第一种情况的信息比第二种情况的信息更有价值。

注:信息量和信息熵的单位都是比特(bit)。

注:在计算信息量或信息熵时,取10的对数lg,或自然常数e的对数ln都是可以的,但是在一次应用过程中,所有的信息量或信息熵都必须采用同一个底。


多维向量空间

一般来说,向量的每个维度之间是不相关的。应尽可能保证维度设置的“正交性”。

例如向量定义:

(姓名,姓,名,出生日期,年龄)

在本例中,“姓名”这个维度可以由“姓”和“名”这两个维度推出,“年龄”也可以由“出生日期”推出。所以说,这种记录方式存在冗余信息,其中一个字段发生变化时,与其相关的其他字段也需要做出变化,这对于保持数据一致性来说,维护成本显然会提高。

在具体场景中,冗余字段也有优点。

例如向量定义:

(用户ID,第一季度销费额,第二季度销费额,第三季度销费额,第四季度销费额, 全年消费总额)

在这种情况下,如果没有“全年消费总额”这一字段,在统计所有用户一年的消费总额时需要将所有的值加起来,在业务反馈时增加了额外的计算量。


[1] 取自《通信的数学理论》,香农,1948。

[2] 哈特莱(R.V.L.Hartley),1928。 


想了解更多关于大数据和机器学习:

分类
程序和算法

归并排序及其Java并发实现

冒泡排序

因为在归并排序的简化过程中需要用到冒泡排序,所以这里先做一下简单介绍。

冒泡排序(Bubble Sort),是一种计算机科学领域的较简单的排序算法。

它重复地走访过要排序的数列,一次比较两个元素,如果他们的顺序错误就把他们交换过来。走访数列的工作是重复地进行直到没有再需要交换,也就是说该数列已经排序完成。

这个算法的名字由来是因为越大的元素会经由交换慢慢“浮”到数列的顶端,故名。

算法原理

  • 比较相邻的元素。如果第一个比第二个大,就交换他们两个。
  • 对每一对相邻元素作同样的工作,从开始第一对到结尾的最后一对。在这一点,最后的元素应该会是最大的数。
  • 针对所有的元素重复以上的步骤,除了最后一个。
  • 持续每次对越来越少的元素重复上面的步骤,直到没有任何一对数字需要比较。

时间复杂度

冒泡排序的最优时间复杂度为:\(O(n)\)

最坏时间复杂度为:\(O(n^2)\)

平均时间复杂度为:\(O(n^2)\)

示例代码

public static void BubbleSort(Integer[] list) {
    for (int i = 0; i < list.length; i++) {
        for (int j = i + 1; j < list.length; j++) {
            if (list[i] > list[j]) {
                swap(list, i, j);
            }
        }
    }
}

在归并排序中用到的是冒泡排序的修改版,多了两个参数first和last,只对数组中索引从first到last的元素进行冒泡排序,示例代码如下。

public static void BubbleSort(
        Integer[] list,
        int first,
        int last) {
    for (int i = first; i <= last; i++) {
        for (int j = i + 1; j <= last; j++) {
            if (list[i] > list[j]) {
                swap(list, i, j);
            }
        }
    }
}

归并排序

归并排序(MERGE-SORT)是建立在归并操作上的一种有效的排序算法,该算法是采用分治法(Divide and Conquer)的一个非常典型的应用。将已有序的子序列合并,得到完全有序的序列;即先使每个子序列有序,再使子序列段间有序。若将两个有序表合并成一个有序表,称为二路归并

算法原理

归并操作的工作原理如下:

  • 第一步:申请空间,使其大小为两个已经排序序列之和,该空间用来存放合并后的序列
  • 第二步:设定两个指针,最初位置分别为两个已经排序序列的起始位置
  • 第三步:比较两个指针所指向的元素,选择相对小的元素放入到合并空间,并移动指针到下一位置
  • 重复步骤3直到某一指针超出序列尾
  • 将另一序列剩下的所有元素直接复制到合并序列尾

时间复杂度

时间复杂度为\(O(nlog₂n)\) 这是该算法中最好、最坏和平均的时间性能。

空间复杂度为 \(O(n)\)

比较操作的次数介于\((nlogn) / 2\)和\(nlogn – n + 1\)。

赋值操作的次数是\((2nlogn)\)。

归并算法的空间复杂度为:\(O (n)\)

归并排序比较占用内存,但却是一种效率高且稳定的算法。

示例代码

public static void MergeArray(
        Integer[] list,
        int first,
        int mid,
        int last) {
    int i = first;
    int j = mid + 1;
    int m = mid;
    int n = last;
    ArrayList<Integer> tmp = new ArrayList<>(
            last - first + 1);
    while (i <= m && j <= n) {
        if (list[i] <= list[j]) {
            tmp.add(list[i++]);
        } else {
            tmp.add(list[j++]);
        }
    }
    while (i <= m) {
        tmp.add(list[i++]);
    }
    while (j <= n) {
        tmp.add(list[j++]);
    }
    for (int index = 0; index < tmp.size(); index++) {
        list[first + index] = tmp.get(index);
    }
}

public static void MergeSort(
    Integer[] list,
    int first,
    int last) {
    if (first < last) {
        int mid = (first + last) / 2;
        MergeSort(list, first, mid);
        MergeSort(list, mid + 1, last);
        MergeArray(list, first, mid, last);
    }
}

归并排序升级版

下图是对1到100个随机数字进行排序,冒泡排序和归并排序运行的时间比较。(单位:秒)

可以看出,虽然归并排序的运行时间小于冒泡排序,但是在排序数字数目较少时(20个以下),冒泡排序的性能要优于归并排序,由此得出以下归并排序的升级版:

设置一个阈值min,当需要排序的数字数目大于min时,继续归并排序;反之,当需要排序的数字数目小于min时,进行冒泡排序。

示例代码如下,这里设置min为需要传入的参数。

public static void MergeSort2(
        Integer[] list,
        int first,
        int last,
        int min) {
    if (last - first > min) {
        int mid = (first + last) / 2;
        MergeSort(list, first, mid);
        MergeSort(list, mid + 1, last);
        MergeArray(list, first, mid, last);
    } else {
        BubbleSort(list, first, last);
    }
}

下面分别对传入参数min为3,4,5,6时对1到100000个随机数字进行排序的运行时间比较。(单位:秒)

放大图像的右下部分。

由于随机数字不稳定的原因,min的不同取值造成的运行时间影响也有所不同。这里的升级版本采取稳定性较好的min=5,示例代码如下。

public static void MergeSort3(
        Integer[] list,
        int first,
        int last) {
    if (last - first > 5) {
        int mid = (first + last) / 2;
        MergeSort(list, first, mid);
        MergeSort(list, mid + 1, last);
        MergeArray(list, first, mid, last);
    } else {
        BubbleSort(list, first, last);
    }
}

归并排序并发实现

对于归并排序考虑到以下问题:

既然把一个排序问题递归地分解成两个更小的排序问题,那么可不可以将这两个排序分别在不同线程中进行呢,答案当然是可以的。

下面使用两种方式进行实现,分别使用Callable/Future和Fork/Join框架,在两个线程中分别进行二路分解,当两个线程都分解完成时,再启用归并程序。(这里都采用上面得到的最优归并排序策略MergeSort3

Callable/Future实现示例代码如下。

public static class MergeSort4 implements Callable<Boolean> {
    private ExecutorService threadPool;
    private Future<Boolean> loadResult0;
    private Future<Boolean> loadResult1;
    private Integer[] list;
    private int first;
    private int last;

    public MergeSort4(Integer[] list, ExecutorService threadPool) {
        this.list = list;
        this.threadPool = threadPool;
        first = 0;
        last = list.length - 1;
    }

    public MergeSort4(Integer[] list, int first, int last,
                      ExecutorService threadPool) {
        this.list = list;
        this.first = first;
        this.last = last;
        this.threadPool = threadPool;
    }

    /**
     * Computes a result, or throws an exception if unable to do so.
     *
     * @return computed result
     * @throws Exception if unable to compute a result
     */
    @Override
    public Boolean call() throws Exception {
        if (last - first > 5) {
            int mid = (first + last) / 2;
            loadResult0 = threadPool.submit(
                    new MergeSort4(list, first, mid, threadPool)
            );
            loadResult1 = threadPool.submit(
                    new MergeSort4(list, mid + 1, last, threadPool)
            );
            try {
                if (loadResult0.get() && loadResult1.get()) {
                    MergeArray(list, first, mid, last);
                }
            } catch (InterruptedException e) {
                e.printStackTrace();
            } catch (ExecutionException e) {
                e.printStackTrace();
            }
        } else {
            BubbleSort(list, first, last);
        }
        return true;
    }
}

Fork/Join实现示例代码如下。

public static class MergeSort5 extends RecursiveTask {
    private Integer[] list;
    private int first;
    private int last;

    public MergeSort5(Integer[] list) {
        this.list = list;
        first = 0;
        last = list.length - 1;
    }

    public MergeSort5(Integer[] list, int first, int last) {
        this.list = list;
        this.first = first;
        this.last = last;
    }

    /**
     * The main computation performed by this task.
     *
     * @return the result of the computation
     */
    @Override
    protected Object compute() {
        if (last - first > 5) {
            int mid = (first + last) / 2;
            MergeSort5 leftTask = new MergeSort5(list, first, mid);
            MergeSort5 rightTask = new MergeSort5(list, mid + 1, last);
            leftTask.fork();
            rightTask.fork();
            leftTask.join();
            rightTask.join();
            MergeArray(list, first, mid, last);
        } else {
            BubbleSort(list, first, last);
        }
        return null;
    }
}

Fork/Join框架详解

Fork/Join框架介绍

Fork/Join框架是Java7提供了的一个用于并行执行任务的框架, 是一个把大任务分割成若干个小任务,最终汇总每个小任务结果后得到大任务结果的框架。

我们再通过Fork和Join这两个单词来理解下Fork/Join框架,Fork就是把一个大任务切分为若干子任务并行的执行,Join就是合并这些子任务的执行结果,最后得到这个大任务的结果。比如计算1+2+……+10000,可以分割成10个子任务,每个子任务分别对1000个数进行求和,最终汇总这10个子任务的结果。当然,也可以分成更多更小的任务。

Fork/Join框架设计

第一步分割任务。首先我们需要有一个fork类来把大任务分割成子任务,有可能子任务还是很大,所以还需要不停的分割,直到分割出的子任务足够小。

第二步执行任务并合并结果。分割的子任务分别放在双端队列里,然后几个启动线程分别从双端队列里获取任务执行。子任务执行完的结果都统一放在一个队列里,启动一个线程从队列里拿数据,然后合并这些数据。

Fork/Join使用两个类来完成以上两件事情:

  • ForkJoinTask:我们要使用ForkJoin框架,必须首先创建一个ForkJoin任务。它提供在任务中执行fork()和join()操作的机制,通常情况下我们不需要直接继承ForkJoinTask类,而只需要继承它的子类,Fork/Join框架提供了以下两个子类:
    • RecursiveAction:用于没有返回结果的任务。
    • RecursiveTask:用于有返回结果的任务。
  • ForkJoinPool:ForkJoinTask需要通过ForkJoinPool来执行,任务分割出的子任务会添加到当前工作线程所维护的双端队列中,进入队列的头部。当一个工作线程的队列里暂时没有任务时,它会随机从其他工作线程的队列的尾部获取一个任务。

Fork/Join框架使用

让我们通过一个简单的需求来使用下Fork/Join框架,需求是:计算1+2+3+4的结果。

使用Fork/Join框架首先要考虑到的是如何分割任务,如果我们希望每个子任务最多执行两个数的相加,那么我们设置分割的阈值是2,由于是4个数字相加,所以Fork/Join框架会把这个任务fork成两个子任务,子任务一负责计算1+2,子任务二负责计算3+4,然后再join两个子任务的结果。

因为是有结果的任务,所以必须继承RecursiveTask,实现代码如下:

public static class CountTask extends RecursiveTask {
    private static final int THRESHOLD = 2;
    private int start;
    private int end;

    public CountTask(int start, int end) {
        this.start = start;
        this.end = end;
    }

    /**
     * The main computation performed by this task.
     *
     * @return the result of the computation
     */
    @Override
    protected Integer compute() {
        int sum = 0;
        boolean canCompute = (end - start) <= THRESHOLD;
        if (canCompute) {
            for (int i = start; i <= end; i++) {
                sum += i;
            }
        } else {
            int middle = (start + end) / 2;
            CountTask leftTask = new CountTask(start, middle);
            CountTask rightTask = new CountTask(middle + 1, end);
            leftTask.fork();
            rightTask.fork();
            int leftResult = (int) leftTask.join();
            int rightResult = (int) rightTask.join();
            sum = leftResult + rightResult;
        }
        return sum;
    }
}

这个类是通用的计算1+2+……+n,每次取两个数相加,当传入参数1和4时,为计算1+2+3+4的结果。


附录1:产生一定数量随机数字的方法

public static Integer[] RandomList(int num) {
    ArrayList<Integer> list = new ArrayList<>(num);
    for (int i = 0; i < num; i++) {
        list.add(i);
    }
    Collections.shuffle(list);
    return list.toArray(new Integer[num]);
}

附录2:交换数组中两个索引的元素的方法

public static void swap(Integer[] list, int index1, int index2) {
    int tmp = list[index1];
    list[index1] = list[index2];
    list[index2] = tmp;
}

附录3:主函数(测试程序)

其中range(delta,maxNum+delta,delta)是测试的数组元素数目的取值集合。

public static void main(String[] args) {
    File file = new File("MergeSort.txt");
    try {
        System.setOut(new PrintStream(file));
    } catch (FileNotFoundException e) {
        e.printStackTrace();
    }

    int maxNum = 100000;
    int delta = 100;
    int listNum = 10000;
    long lastTime, curTime, nsPerTime;
    Integer[] list;

    System.out.println("BubbleSort");
    for (listNum = delta; listNum <= maxNum; listNum += delta) {
        list = RandomList(listNum);
        lastTime = System.nanoTime();
        BubbleSort(list);
        curTime = System.nanoTime();
        nsPerTime = curTime - lastTime;
        System.out.print(nsPerTime / 1.0E9 + ",");
    }
    System.out.println();

    System.out.println("MergeSort");
    for (listNum = delta; listNum <= maxNum; listNum += delta) {
        list = RandomList(listNum);
        lastTime = System.nanoTime();
        MergeSort(list, 0, list.length - 1);
        curTime = System.nanoTime();
        nsPerTime = curTime - lastTime;
        System.out.print(nsPerTime / 1.0E9 + ",");
    }
    System.out.println();

    System.out.println("MergeSort 3");
    for (listNum = delta; listNum <= maxNum; listNum += delta) {
        list = RandomList(listNum);
        lastTime = System.nanoTime();
        MergeSort2(list, 0, list.length - 1, 3);
        curTime = System.nanoTime();
        nsPerTime = curTime - lastTime;
        System.out.print(nsPerTime / 1.0E9 + ",");
    }
    System.out.println();

    System.out.println("MergeSort 4");
    for (listNum = delta; listNum <= maxNum; listNum += delta) {
        list = RandomList(listNum);
        lastTime = System.nanoTime();
        MergeSort2(list, 0, list.length - 1, 4);
        curTime = System.nanoTime();
        nsPerTime = curTime - lastTime;
        System.out.print(nsPerTime / 1.0E9 + ",");
    }
    System.out.println();

    System.out.println("MergeSort 5");
    for (listNum = delta; listNum <= maxNum; listNum += delta) {
        list = RandomList(listNum);
        lastTime = System.nanoTime();
        MergeSort2(list, 0, list.length - 1, 5);
        curTime = System.nanoTime();
        nsPerTime = curTime - lastTime;
        System.out.print(nsPerTime / 1.0E9 + ",");
    }
    System.out.println();

    System.out.println("MergeSort 6");
    for (listNum = delta; listNum <= maxNum; listNum += delta) {
        list = RandomList(listNum);
        lastTime = System.nanoTime();
        MergeSort2(list, 0, list.length - 1, 6);
        curTime = System.nanoTime();
        nsPerTime = curTime - lastTime;
        System.out.print(nsPerTime / 1.0E9 + ",");
    }
    System.out.println();

    System.out.println("MergeSort Callable");
    for (listNum = delta; listNum <= maxNum; listNum += delta) {
        list = RandomList(listNum);
        ExecutorService threadPool = Executors.newCachedThreadPool();
        MergeSort4 mergeSort4 = new MergeSort4(list, threadPool);
        lastTime = System.nanoTime();
        Future<Boolean> result = threadPool.submit(mergeSort4);
        try {
            result.get();
        } catch (InterruptedException e) {
            e.printStackTrace();
        } catch (ExecutionException e) {
            e.printStackTrace();
        } finally {
            threadPool.shutdown();
        }
        curTime = System.nanoTime();
        nsPerTime = curTime - lastTime;
        System.out.print(nsPerTime + ",");
    }
    System.out.println();

    System.out.println("MergeSort RecursiveTask");
    for (listNum = delta; listNum <= maxNum; listNum += delta) {
        list = RandomList(listNum);
        ForkJoinPool forkJoinPool = new ForkJoinPool();
        MergeSort5 mergeSort5 = new MergeSort5(list);
        lastTime = System.nanoTime();
        Future result = forkJoinPool.submit(mergeSort5);
        try {
            result.get();
        } catch (InterruptedException e) {
            e.printStackTrace();
        } catch (ExecutionException e) {
            e.printStackTrace();
        } finally {
            forkJoinPool.shutdown();
        }
        curTime = System.nanoTime();
        nsPerTime = curTime - lastTime;
        System.out.print(nsPerTime + ",");
    }
    System.out.println();
}

附录4:源代码

import java.io.File;
import java.io.FileNotFoundException;
import java.io.PrintStream;
import java.util.ArrayList;
import java.util.Collections;
import java.util.concurrent.*;

public class MergeSort {
    public static void main(String[] args) {
        File file = new File("MergeSort.txt");
        try {
            System.setOut(new PrintStream(file));
        } catch (FileNotFoundException e) {
            e.printStackTrace();
        }

        int maxNum = 100000;
        int delta = 100;
        int listNum = 10000;
        long lastTime, curTime, nsPerTime;
        Integer[] list;

        System.out.println("BubbleSort");
        for (listNum = delta; listNum <= maxNum; listNum += delta) {
            list = RandomList(listNum);
            lastTime = System.nanoTime();
            BubbleSort(list);
            curTime = System.nanoTime();
            nsPerTime = curTime - lastTime;
            System.out.print(nsPerTime / 1.0E9 + ",");
        }
        System.out.println();

        System.out.println("MergeSort");
        for (listNum = delta; listNum <= maxNum; listNum += delta) {
            list = RandomList(listNum);
            lastTime = System.nanoTime();
            MergeSort(list, 0, list.length - 1);
            curTime = System.nanoTime();
            nsPerTime = curTime - lastTime;
            System.out.print(nsPerTime / 1.0E9 + ",");
        }
        System.out.println();

        System.out.println("MergeSort 3");
        for (listNum = delta; listNum <= maxNum; listNum += delta) {
            list = RandomList(listNum);
            lastTime = System.nanoTime();
            MergeSort2(list, 0, list.length - 1, 3);
            curTime = System.nanoTime();
            nsPerTime = curTime - lastTime;
            System.out.print(nsPerTime / 1.0E9 + ",");
        }
        System.out.println();

        System.out.println("MergeSort 4");
        for (listNum = delta; listNum <= maxNum; listNum += delta) {
            list = RandomList(listNum);
            lastTime = System.nanoTime();
            MergeSort2(list, 0, list.length - 1, 4);
            curTime = System.nanoTime();
            nsPerTime = curTime - lastTime;
            System.out.print(nsPerTime / 1.0E9 + ",");
        }
        System.out.println();

        System.out.println("MergeSort 5");
        for (listNum = delta; listNum <= maxNum; listNum += delta) {
            list = RandomList(listNum);
            lastTime = System.nanoTime();
            MergeSort2(list, 0, list.length - 1, 5);
            curTime = System.nanoTime();
            nsPerTime = curTime - lastTime;
            System.out.print(nsPerTime / 1.0E9 + ",");
        }
        System.out.println();

        System.out.println("MergeSort 6");
        for (listNum = delta; listNum <= maxNum; listNum += delta) {
            list = RandomList(listNum);
            lastTime = System.nanoTime();
            MergeSort2(list, 0, list.length - 1, 6);
            curTime = System.nanoTime();
            nsPerTime = curTime - lastTime;
            System.out.print(nsPerTime / 1.0E9 + ",");
        }
        System.out.println();

        System.out.println("MergeSort Callable");
        for (listNum = delta; listNum <= maxNum; listNum += delta) {
            list = RandomList(listNum);
            ExecutorService threadPool = Executors.newCachedThreadPool();
            MergeSort4 mergeSort4 = new MergeSort4(list, threadPool);
            lastTime = System.nanoTime();
            Future<Boolean> result = threadPool.submit(mergeSort4);
            try {
                result.get();
            } catch (InterruptedException e) {
                e.printStackTrace();
            } catch (ExecutionException e) {
                e.printStackTrace();
            } finally {
                threadPool.shutdown();
            }
            curTime = System.nanoTime();
            nsPerTime = curTime - lastTime;
            System.out.print(nsPerTime + ",");
        }
        System.out.println();

        System.out.println("MergeSort RecursiveTask");
        for (listNum = delta; listNum <= maxNum; listNum += delta) {
            list = RandomList(listNum);
            ForkJoinPool forkJoinPool = new ForkJoinPool();
            MergeSort5 mergeSort5 = new MergeSort5(list);
            lastTime = System.nanoTime();
            Future result = forkJoinPool.submit(mergeSort5);
            try {
                result.get();
            } catch (InterruptedException e) {
                e.printStackTrace();
            } catch (ExecutionException e) {
                e.printStackTrace();
            } finally {
                forkJoinPool.shutdown();
            }
            curTime = System.nanoTime();
            nsPerTime = curTime - lastTime;
            System.out.print(nsPerTime + ",");
        }
        System.out.println();
    }

    public static Integer[] combine(Integer[] a, Integer[] b) {
        ArrayList<Integer> list = new ArrayList<>(
                a.length + b.length);
        int i = 0;
        int j = 0;
        while (i < a.length && j < b.length) {
            if (a[i] < b[j]) {
                list.add(a[i++]);
            } else {
                list.add(b[j++]);
            }
        }
        while (i < a.length) {
            list.add(a[i++]);
        }
        while (j < b.length) {
            list.add(b[j++]);
        }
        return list.toArray(new Integer[a.length + b.length]);
    }

    public static void MergeArray(
            Integer[] list,
            int first,
            int mid,
            int last) {
        int i = first;
        int j = mid + 1;
        int m = mid;
        int n = last;
        ArrayList<Integer> tmp = new ArrayList<>(
                last - first + 1);
        while (i <= m && j <= n) {
            if (list[i] <= list[j]) {
                tmp.add(list[i++]);
            } else {
                tmp.add(list[j++]);
            }
        }
        while (i <= m) {
            tmp.add(list[i++]);
        }
        while (j <= n) {
            tmp.add(list[j++]);
        }
        for (int index = 0; index < tmp.size(); index++) {
            list[first + index] = tmp.get(index);
        }
    }

    public static void MergeSort(
            Integer[] list,
            int first,
            int last) {
        if (first < last) {
            int mid = (first + last) / 2;
            MergeSort(list, first, mid);
            MergeSort(list, mid + 1, last);
            MergeArray(list, first, mid, last);
        }
    }

    public static void MergeSort2(
            Integer[] list,
            int first,
            int last,
            int min) {
        if (last - first > min) {
            int mid = (first + last) / 2;
            MergeSort(list, first, mid);
            MergeSort(list, mid + 1, last);
            MergeArray(list, first, mid, last);
        } else {
            BubbleSort(list, first, last);
        }
    }

    public static void MergeSort3(
            Integer[] list,
            int first,
            int last) {
        if (last - first > 5) {
            int mid = (first + last) / 2;
            MergeSort(list, first, mid);
            MergeSort(list, mid + 1, last);
            MergeArray(list, first, mid, last);
        } else {
            BubbleSort(list, first, last);
        }
    }

    public static class MergeSort4 implements Callable<Boolean> {
        private ExecutorService threadPool;
        private Future<Boolean> loadResult0;
        private Future<Boolean> loadResult1;
        private Integer[] list;
        private int first;
        private int last;
        public MergeSort4(Integer[] list, ExecutorService threadPool) {
            this.list = list;
            this.threadPool = threadPool;
            first = 0;
            last = list.length - 1;
        }
        public MergeSort4(Integer[] list, int first, int last,
                          ExecutorService threadPool) {
            this.list = list;
            this.first = first;
            this.last = last;
            this.threadPool = threadPool;
        }

        /**
         * Computes a result, or throws an exception if unable to do so.
         *
         * @return computed result
         * @throws Exception if unable to compute a result
         */
        @Override
        public Boolean call() throws Exception {
            if (last - first > 5) {
                int mid = (first + last) / 2;
                loadResult0 = threadPool.submit(
                        new MergeSort4(list, first, mid, threadPool)
                );
                loadResult1 = threadPool.submit(
                        new MergeSort4(list, mid + 1, last, threadPool)
                );
                try {
                    if (loadResult0.get() && loadResult1.get()) {
                        MergeArray(list, first, mid, last);
                    }
                } catch (InterruptedException e) {
                    e.printStackTrace();
                } catch (ExecutionException e) {
                    e.printStackTrace();
                }
            } else {
                BubbleSort(list, first, last);
            }
            return true;
        }
    }

    public static class MergeSort5 extends RecursiveTask {
        private Integer[] list;
        private int first;
        private int last;

        public MergeSort5(Integer[] list) {
            this.list = list;
            first = 0;
            last = list.length - 1;
        }

        public MergeSort5(Integer[] list, int first, int last) {
            this.list = list;
            this.first = first;
            this.last = last;
        }

        /**
         * The main computation performed by this task.
         *
         * @return the result of the computation
         */
        @Override
        protected Object compute() {
            if (last - first > 5) {
                int mid = (first + last) / 2;
                MergeSort5 leftTask = new MergeSort5(list, first, mid);
                MergeSort5 rightTask = new MergeSort5(list, mid + 1, last);
                leftTask.fork();
                rightTask.fork();
                leftTask.join();
                rightTask.join();
                MergeArray(list, first, mid, last);
            } else {
                BubbleSort(list, first, last);
            }
            return null;
        }
    }

    public static class CountTask extends RecursiveTask {
        private static final int THRESHOLD = 2;
        private int start;
        private int end;

        public CountTask(int start, int end) {
            this.start = start;
            this.end = end;
        }

        /**
         * The main computation performed by this task.
         *
         * @return the result of the computation
         */
        @Override
        protected Integer compute() {
            int sum = 0;
            boolean canCompute = (end - start) <= THRESHOLD;
            if (canCompute) {
                for (int i = start; i <= end; i++) {
                    sum += i;
                }
            } else {
                int middle = (start + end) / 2;
                CountTask leftTask = new CountTask(start, middle);
                CountTask rightTask = new CountTask(middle + 1, end);
                leftTask.fork();
                rightTask.fork();
                int leftResult = (int) leftTask.join();
                int rightResult = (int) rightTask.join();
                sum = leftResult + rightResult;
            }
            return sum;
        }
    }

    public static Integer[] RandomList(int num) {
        ArrayList<Integer> list = new ArrayList<>(num);
        for (int i = 0; i < num; i++) {
            list.add(i);
        }
        Collections.shuffle(list);
        return list.toArray(new Integer[num]);
    }

    public static void BubbleSort(Integer[] list) {
        for (int i = 0; i < list.length; i++) {
            for (int j = i + 1; j < list.length; j++) {
                if (list[i] > list[j]) {
                    swap(list, i, j);
                }
            }
        }
    }

    public static void BubbleSort(
            Integer[] list,
            int first,
            int last) {
        for (int i = first; i <= last; i++) {
            for (int j = i + 1; j <= last; j++) {
                if (list[i] > list[j]) {
                    swap(list, i, j);
                }
            }
        }
    }

    public static void swap(Integer[] list, int index1, int index2) {
        int tmp = list[index1];
        list[index1] = list[index2];
        list[index2] = tmp;
    }
}