• Archive by category "R语言实践"

Blog Archives

R语言实现聚类kmeans

R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。

要成为有理想的极客,我们不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让我们一起动起来吧,开始R的极客理想。

关于作者:

  • 张丹(Conan), 程序员/Quant: Java,R,Nodejs
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-cluster-kmeans

前言

聚类属于无监督学习中的一种方法,k-means作为数据挖掘的十大算法之一,是一种最广泛使用的聚类算法。我们使用聚类算法将数据集的点,分到特定的组中,同一组的数据点具有相似的特征,而不同类中的数据点特征差异很大。PAM是对k-means的一种改进算法,能降低异常值对于聚类效果的影响。

聚类可以帮助我们认识未知的数据,发现新的规律。

目录

  1. k-means实现
  2. PAM实现
  3. 可视化和段剖面图

1. k-means实现

k-means算法,是一种最广泛使用的聚类算法。k-means以k作为参数,把数据分为k个组,通过迭代计算过程,将各个分组内的所有数据样本的均值作为该类的中心点,使得组内数据具有较高的相似度,而组间的相似度最低。

k-means工作原理:

  1. 初始化数据,选择k个对象作为中心点。
  2. 遍历整个数据集,计算每个点与每个中心点的距离,将它分配给距离中心最近的组。
  3. 重新计算每个组的平均值,作为新的聚类中心。
  4. 上面2-3步,过程不断重复,直到函数收敛,不再新的分组情况出现。

k-means聚类,适用于连续型数据集。在计算数据样本之间的距离时,通常使用欧式距离作为相似性度量。k-means支持多种距离计算,还包括maximum, manhattan, pearson, correlation, spearman, kendall等。各种的距离算法的介绍,请参考文章R语言实现46种距离算法

1.1 kmeans()函数实现

在R语言中,我们可以直接调用系统中自带的kmeans()函数,就可以实现k-means的聚类。同时,有很多第三方算法包也提供了k-means的计算函数。当我们需要使用kmeans算法,可以使用第三方扩展的包,比如flexclust, amap等包。

本文的系统环境为:

  • Win10 64bit
  • R: 3.4.4 x86_64-w64-mingw32

接下来,让我们做一个k-means聚类的例子。首先,创建数据集。

# 创建数据集
> set.seed(0)
> df <- rbind(matrix(rnorm(100, 0.5, 4.5), ncol = 2),
+             matrix(rnorm(100, 0.5, 0.1), ncol = 2))
> colnames(df) <- c("x", "y")
> head(df)
              x          y
[1,]  6.1832943  1.6976181
[2,] -0.9680501 -1.1951622
[3,]  6.4840967 11.4861408
[4,]  6.2259319 -3.0790260
[5,]  2.3658865  0.2530514
[6,] -6.4297752  1.6256360

使用stats::kmeans()函数,进行聚类。


> cl <- kmeans(df,2); cl
K-means clustering with 2 clusters of sizes 14, 86

Cluster means:                   # 中心点坐标
          x         y
1  5.821526 2.7343127
2 -0.315946 0.1992429

Clustering vector:               # 分组的索引
  [1] 1 2 1 1 2 2 2 2 2 1 1 2 2 2 2 2 2 2 2 2 2 2 2 1 2 1 1 2 2 2 2 2 2 2 2 1 1 2 1 2 1 2 2 2 2 2 2 1 1 2
 [51] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

Within cluster sum of squares by cluster:   
[1] 316.0216 716.4009                       # withinss,分组内平方和  
 (between_SS / total_SS =  34.0 %)          # 组间的平方和/总平方和,用于衡量点聚集程度

Available components:            # 对象属性
[1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss" "betweenss"   
[7] "size"         "iter"         "ifault"      

# 查看数据分组情况,第1组86个,第2组14个
> cl$size
[1] 86 14

对象属性解读:

  • cluster,每个点的分组
  • centers,聚类的中心点坐标
  • totss,总平方和
  • withinss,每个分组内的平方和
  • tot.withinss,分组总和,sum(withinss)
  • betweenss,组间的平方和,totss – tot.withinss
  • size,每个组中的数据点数量
  • iter,迭代次数。
  • ifault,可能有问题的指标

1.2 kcca()函数实现
我们再使用flexclust::kcca()函数,进行聚类。


# 安装flexclust包
> # install.packages("flexclust")
> library(flexclust)

# 进行聚类
> clk<-kcca(df,k=2);clk
kcca object of family ‘kmeans’ 

call:
kcca(x = df, k = 2)

cluster sizes:  # 聚类的分组大小
 1  2 
84 16 

# 聚类的中心
> clk@centers
              x         y
[1,] -0.3976465 0.2015319
[2,]  5.4832702 2.4054118

# 查看聚类的概览信息
> summary(clk)
kcca object of family ‘kmeans’ 

call:
kcca(x = df, k = 2)

cluster info:         # 每个组的基本信息,包括分组数量,平均距离、最大距离、分割值
  size  av_dist max_dist separation
1   84 2.102458 9.748136   3.368939
2   16 3.972920 9.576635   3.189891

convergence after 5 iterations                   # 5次迭代
sum of within cluster distances: 240.1732        # 聚类距离之和

我们比较2个不同包的k-means算法,所得到的分组数据都是一样的,中心点位置略有一点偏差。接下来,我们可以把聚类画图。

> plot(df, col = cl$cluster, main="Kmeans Cluster")
> points(cl$centers, col = 1:3, pch = 10, cex = 4) # 画出kmeans()函数效果

从上图中看到k-means的总分2组,每个组的中心点分别用红色十字圆圈和黑色十字圆圈表示,为组内的所有数据样本的均值。再叠加上kcca()函数聚类后的中心点画图。

> points(clk@centers, col = 3:4, pch = 10, cex = 4)  # 画出kcca()函数效果


新的中心点,分别用别用绿色十字圆圈和蓝色十字圆圈表示。虽然我们使用了相同的算法,分组个数也相同,但中心点还有一些不同的。

这里其实就要对聚类的稳定性进行判断了,有可能是聚类迭代次数过少,就会出现不同的聚类结果,就需要增加迭代次数,达到每次计算结果是一致的。也有可能是因为不同的包,实现的代码有所区别导致的。

k-means算法,也有一些缺点就是对于孤立点是敏感的,会被一些极端值影响聚类的效果。一种改进的算法是PAM,用于解决这个问题。PAM不使用分组平均值作为计算的参照点,而是直接使用每个组内最中心的对象作为中心点。

2. PAM实现

PAM(Partitioning Around Medoids),又叫k-medoids,它可以将数据分组为k个组,k为数量是要事前定义的。PAM与k-means一样,找到距离中心点最小点组成同一类。PAM对噪声和异常值更具鲁棒性,该算法的目标是最小化对象与其最接近的所选对象的平均差异。PAM可以支持混合的数据类型,不仅限于连续变量。

PAM算法分为两个阶段:

  1. 第1阶段BUILD,为初始集合S选择k个对象的集合。
  2. 第2阶段SWAP,尝试用未选择的对象,交换选定的中心点,来提高聚类的质量。

PAM的工作原理:

  1. 初始化数据集,选择k个对象作为中心。
  2. 遍历数据点,把每个数据点关联到最近中心点m。
  3. 随机选择一个非中心对象,与中心对象交换,计算交换后的距离成本
  4. 如果总成本增加,则撤销交换的动作。
  5. 上面2-4步,过程不断重复,直到函数收敛,中心不再改变为止。

优点与缺点:

  • 消除了k-means算法对于孤立点的敏感性。
  • 比k-means的计算的复杂度要高。
  • 与k-means一样,必须设置k的值。
  • 对小的数据集非常有效,对大数据集效率不高。

在R语言中,我们可以通过cluster包来使用pam算法函数。cluster包的安装很简单,一条命令就安装完了。


> install.packages("cluster")
> library(cluster)

pam()函数定义:


pam(x, k, diss = inherits(x, "dist"), metric = "euclidean",
    medoids = NULL, stand = FALSE, cluster.only = FALSE,
    do.swap = TRUE,
    keep.diss = !diss && !cluster.only && n < 100,
    keep.data = !diss && !cluster.only,
    pamonce = FALSE, trace.lev = 0)

参数列表:

  • x,数据框或矩阵,允许有空值(NA)
  • k,设置分组数量
  • diss,为TRUE时,x为距离矩阵;为FALSE时,x是变量矩阵。默认为FALSE
  • metric,设置距离算法,默认为euclidean,距离矩阵忽略此项
  • medoids,指定初始的中心,默认为不指定。
  • stand,为TRUE时进行标准化,距离矩阵忽略此项。
  • cluster.only,为TRUE时,仅计算聚类结果,默认为FALSE
  • do.swap,是否进行中心点交换,默认为TRUE;对于超大的数据集,可以不进行交换。
  • keep.diss,是否保存距离矩阵数据
  • keep.data,是否保存原始数据
  • pamonce,一种加速算法,接受值为TRUE,FALSE,0,1,2
  • trace.lev,日志打印,默认为0,不打印

我们使用上面已创建好的数据集df,进行pam聚类,设置k=2。

> kclus <- pam(df,2)

# 查看kclus对象
> kclus
Medoids:                                     # 中心点
     ID         x         y
[1,] 27 5.3859621 1.1469717
[2,] 89 0.4130217 0.4798659

Clustering vector:                           # 分组
  [1] 1 2 1 1 2 2 2 2 2 1 1 2 2 2 2 2 2 2 2 2 2 2 2 1 2 1 1 2 2 2 2 2 2 2 1 1 1 2 1 2 1 2 2 2 2 2 2 1 1 2
 [51] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

Objective function:                          # 目标函数的局部最小值
   build     swap                           
2.126918 2.124185 

Available components:                        # 聚类对象的属性
 [1] "medoids"    "id.med"     "clustering" "objective"  "isolation"  "clusinfo"   "silinfo"   
 [8] "diss"       "call"       "data"      

> kclus$clusinfo        # 聚类的分组数量,每个组的平均距离、最大距离、分割值
     size  max_diss  av_diss diameter separation
[1,]   15 10.397323 4.033095 17.35984   1.556862
[2,]   85  9.987604 1.787318 15.83646   1.556862

属性解读:

  • medoids,中心点的数据值
  • id.med,中心点的索引
  • clustering,每个点的分组
  • objective,目标函数的局部最小值
  • isolation,孤立的聚类(用L或L*表示)
  • clusinfo,每个组的基本信息
  • silinfo,存储各观测所属的类、其邻居类以及轮宽(silhouette)值
  • diss,不相似度
  • call,执行函数和参数
  • data,原始数据集

把聚类画图输出。

# 画图
> plot(df, col = kclus$clustering, main="Kmedoids Cluster")
> points(kclus$medoids, col = 1:3, pch = 10, cex = 4)

图中,PAM聚类后分为2组,红色一组,黑色一组,用十字圆圈表示2个中心点,可以清晰地看到中心点就是数据点。

我们可以在开始计算时,设置聚类的中心点,为索引1,2坐标点,打印聚类的日志,查看计算过程。


# 设置聚类的中心为1,2
> kclus2<-pam(df,2,medoids=c(1,2),trace.lev=20)
C pam(): computing 4951 dissimilarities from  100 x 2  matrix: [Ok]
pam()'s bswap(*, s=21.837, pamonce=0): medoids given
  after build: medoids are   1   2
  and min.dist dysma[1:n] are
      0      0   9.79   4.78   3.63   6.15   5.23  0.929   8.44   8.59
   2.29   2.69   4.48   1.19   1.98   2.81   5.39    4.2   3.72   4.56
   1.84   3.99    2.4    2.7   4.84   5.08  0.969   2.01   4.94   5.06
   1.94    7.4   5.19   1.62   3.94   3.12   3.51   0.65   4.46   4.61
   5.16   4.57   1.82   3.21   5.79   4.01   5.59   5.38   1.95    6.2
   2.41   2.09    2.2   2.43   2.24   2.26   2.09   2.39   2.21   2.33
   2.24   2.14   2.45   2.37    2.2   2.37   2.13   2.33   2.25   2.18
   2.38   2.19   2.15   2.14    2.1   2.39   2.24   2.24   2.12   2.14
   2.34   2.18   2.25   2.26   2.33   2.17   2.18   2.12   2.17   2.27
   2.29   2.26   2.38   2.12   2.25   2.33   2.09   2.21   2.24   2.13
   swp new  89 <->   2 old; decreasing diss. 306.742 by -93.214
   swp new  27 <->   1 old; decreasing diss. 213.528 by -1.10916
end{bswap()}, end{cstat()}

# 查看中心
> kclus2$id.med
[1] 27 89

通过日志查看,我们可以清楚地看到,2个中心的选择过程,分别用89替换1,距离成本减少93.214,用27替换2,距离成本减少1.1。

PAM作为k-means的一种改进算法,到底结果是否更合理,还要看最终哪种结果能够准确地表达业务的含义,被业务人员所认可,就需要不断地和业务人员来沟通。

3. 可视化和段剖面图

我们实现了聚类计算后,通常需要把复杂的数据逻辑,用简单的语言和图形来解释给业务人员,聚类的可视化就很重要的。如果数据量不太大,参与聚类的指标维度不太多的时候,我们可以用2维散点图,把指标两两画出来。

我们对iris数据集,进行k-means聚类分成3组,画出聚类后的2维散点图结果。

> res <- kmeans(iris[,1:4], centers=3)
> pairs(iris, col = res$cluster + 1)


每2个维度就会生成一张图, 我们可以全面直观的看到聚类的效果。

高级画图工具,使用GGally包中的ggpairs()函数。

> library(GGally)
> ggpairs(iris,columns = 1:5,mapping=aes(colour=as.character(res$cluster)))


图更漂亮了而且包含更多的信息,除了2维散点图,还包括了相关性检查,分布图,分箱图,频率图等。用这样的可视化效果图与业务人员沟通,一定会非常愉快的。

但是如果数据维度,不止3个而是30个,数据量也不是几百个点,而是几百万个点,再用2维散点图画出来就会很难看了,而且也表达不清,还会失去重点,计算的复杂度也是非常的高。

当数据量和数据维度多起来,我们就需要用段剖面图来做展示了,放弃个体特征,反应的群体特征和规律。

使用flexclust包中的barchart()函数,画出段剖面图,我们还是用iris数据集进行举例。


> library(flexclust)
> clk2 <- cclust(iris[,-5], k=3);clk2
kcca object of family ‘kmeans’ 

call:
cclust(x = iris[, -5], k = 3)

cluster sizes:
 1  2  3 
39 61 50 

# 画出段剖面图
> barchart(clk2,legend=TRUE)

如上图所示,每一区块是一个类别,每行是不同的指标。红点表示均值,柱状是这个类别每个指标的情况,透明色表示不重要指标。

查看段剖面图,可以清楚的看到,每个分组中特征是非常明显的。

  • Cluster1中,有39个数据点占26%,Sepal.Width指标在均值附近,其他指标都大于均值。
  • Cluster2中,有61个数据点占41%,Sepal.Width指标略小于均值,其他指标在均值附近。
  • Cluster3中,有50个数据点占33%,Sepal.Width略大于均值,其他指标都小于均值。

从段剖面图,我们可以一眼就能直观地发现数据聚类后的每个分组的总体特征,而不是每个分组中数据的个体特征,对于数据的解读是非常有帮助的。

对于段剖面图,原来我并不知道是什么效果。在和业务人员沟通中,发现他们使用SAS软件做出了很漂亮的段剖面图,而且他们都能理解,后来我发现R语言也有这个工具函数,图确实能极大地帮助进行数据解读,所以写了这篇文章记录一下。

本文介绍了k-means的聚类计算方法和具体的使用方法,也是对最近做了一个聚类模型的总结。作为数据分析师,我们不仅自己能发现数据的规律,还要让业务人员看明白你的思路,看懂数据的价值,这也是算法本身的价值。

转载请注明出处:
http://blog.fens.me/r-cluster-kmeans

打赏作者

用R语言实现信息度量

R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。

要成为有理想的极客,我们不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让我们一起动起来吧,开始R的极客理想。

关于作者:

  • 张丹(Conan), 程序员/Quant: Java,R,Nodejs
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-entropy

前言

香农的《通信的数学理论》是20世纪非常伟大的著作,被认为是现代信息论研究的开端。信息论定义了信息熵,用于把信息进行度量,以比特(bit)作为量纲单位,为如今发达的信息产业和互联网产业奠定了基础。本文接上一篇文章R语言实现46种距离算法,继续philentropy包的介绍,包括信息度量函数的使用。

目录

  1. 信息熵介绍
  2. 关键概念
  3. 信息度量函数
  4. 应用举例

1.信息熵介绍

信息论(Information Theory)是概率论与数理统计的一个分枝,用于研究信息处理、信息熵、通信系统、数据传输、率失真理论、密码学、信噪比、数据压缩等问题的应用数学学科。信息论将信息的传递作为一种统计现象来考虑,给出了估算通信信道容量的方法。信息传输和信息压缩是信息论研究中的两大领域。

香农被称为是“信息论之父”,香农于1948年10月发表的A Mathematical Theory of Communication,通信的数学理论(中文版),通常被认为是现代信息论研究的开端。

信息熵,是对信息随机性的量度,又指信息能被压缩的极限,用bit作为衡量信息的最小单位。一切信息所包含的信息量,都是1bit的正整数倍。计算机系统中常采用二进制编码,一个0或1就是1bit。

举例来说明一下信息熵的计算原理,假设小明最喜欢5种水果,苹果、香蕉、西瓜、草莓、樱桃中的一种,如果小明没有偏爱,选择每种水果的概率都是20%,那么这一信息的信息熵为

H(A) = -1*(0.2*log2(0.2)*5)
= 2.321928 bits

如果小明偏爱香蕉,选择这5种水果的概率分别是10%,20%,45%,15%,10%,那么这一信息信息熵为


H(B)=-1*(0.1*log2(0.1)+0.2*log2(0.2)+0.45*log2(0.45)+0.15*log2(0.15)+0.1*log2(0.1))
= 2.057717 bits

从结果得到H(A)大于H(B),信息熵越大表示越不确定。对于B的情况,对某一种水果的偏好,比A增加了确定性的因素,所以H(B)小于H(A)是符合对于信息熵的定义的。

2. 关键概念

我们从一幅图来认识信息熵,图中显示了随机变量X和Y的2个集合,在信息熵的概念里的所有可能逻辑关系。两个圆所包含的面积为联合熵H(X,Y), 左边的整个圆表示X的熵H(X),左边半圆是条件熵H(X|Y)。 右边的整个圆表示Y的熵H(Y),右边半圆条件熵H(Y|X),中间交集的部分是互信息I(X; Y)

信息熵(Entropy):是对信息随机性的量度,用于计算信息能被压缩的极限。对随机变量X,不确定性越大,X的信息熵H(X)也就越大。

公式定义:

H(x)的取值范围,0<=H(x)<=log(n), 其中n是随机变量x取值的种类数。需要注意的是,熵只依赖于随机变量的分布,与随机变量取值无关。

条件熵(Conditional Entropy):表示两个随机变量X和Y,在已知Y的情况下对随机变量X的不确定性,称之为条件熵H(X|Y),

公式定义:

联合熵(Joint Entropy):表示为两个随机事件X和Y的熵的并集,联合熵解决将一维随机变量分布推广到多维随机变量分布。

公式定义:

互信息(Mutual Information, 信息增益):两个随机变量X和Y,Y对X的互信息,为后验概率与先验概率比值的对数,即原始的熵H(X)和已知Y的情况下的条件熵H(X|Y)的比值的对数,信息增益越大表示条件Y对于确定性的贡献越大。互信息,也可以用来衡量相似性。

公式定义:

当MI(X,Y)=0时,表示两个事件X和Y完全不相关。决策树ID3算法就是使用信息增益来划分特征,信息增益大时,说明对数据划分帮助很大,优先选择该特征进行决策树的划分。

信息增益比率:是信息增益与该特征的信息熵之比,用于解决信息增益对多维度特征的选择,决策树C4.5算法使用信息增益比率进行特征划分。

KL散度(Kullback–Leibler Divergence, 相对熵):随机变量x取值的两个概率分布p和q,用来衡量这2个分布的差异,通常用p表示真实分布,用q表示预测分布。

公式定义:

n为事件的所有可能性,如果两个分布完全相同,那么它们的相关熵为0。如果相对熵KL越大,说明它们之间的差异越大,反之相对熵KL越小,说明它们之间的差异越小。

交叉熵(Cross Entropy):是对KL散度的一种变型,把KL散度log(p(x)/q(x))进行拆分,前面部分就是p的熵H(p),后面就是交叉熵H(p,q)。

公式定义:

交叉熵可以用来计算学习模型分布与训练分布之间的差异,一般在机器学习中直接用交叉熵做损失函数,用于评估模型。

信息论是通信理论的基础,也是xx的基础,关于信息论的理论,等后面有时时间再做分享,本文重要研究信息熵的函数计算问题。

3. 信息度量函数

philentropy包的函数,主要分为3种类别的函数,第一类是距离测量的函数,第二类是相关性分析,第三类是信息度量函数,本文重点介绍这些信息度量的函数。有关于距离测量函数和相关性分析函数,请参考文章R语言实现46种距离算法

我们来看一下,philentropy包里信息度量的函数:

  • H(): 香农熵, Shannon’s Entropy H(X)
  • JE() : 联合熵, Joint-Entropy H(X,Y)
  • CE() : 条件熵, Conditional-Entropy H(X|Y)
  • MI() : 互信息, Shannon’s Mutual Information I(X,Y)
  • KL() : KL散度, Kullback–Leibler Divergence
  • JSD() : JS散度,Jensen-Shannon Divergence
  • gJSD() : 通用JS散度,Generalized Jensen-Shannon Divergence

本文的系统环境为:

  • Win10 64bit
  • R: 3.4.2 x86_64-w64-mingw32

3.1 H()香农熵
H()函数,可用于快速计算任何给定概率向量的香农熵。

H()函数定义:

H (x, unit = "log2") 

参数列表:

  • x, 概率向量
  • unit,对数化的单位,默认为log2

函数使用:


# 创建数据x
> x<-1:10;x
 [1]  1  2  3  4  5  6  7  8  9 10
> x1<-x/sum(x);x1
 [1] 0.01818182 0.03636364 0.05454545 0.07272727
 [5] 0.09090909 0.10909091 0.12727273 0.14545455
 [9] 0.16363636 0.18181818

# 计算香农熵
> H(px)
[1] 3.103643

同样地,我们也可以用程序实现公式自己算一下。


# 创建数据x
> x<-1:10
#计算x的概率密度px
> px<-x/sum(x)  

# 根据公式计算香农熵
> -1*sum(px*log2(px))
[1] 3.103643

我们动手的计算结果,用于H()函数的计算结果是一致的。

3.2 CE()条件熵

CE()函数,基于给定的联合概率向量P(X,Y)和概率向量P(Y),根据公式 H(X|Y)= H(X,Y)-H(Y)计算香农的条件熵。

函数定义:

CE(xy, y, unit = "log2")

参数列表:

  • xy, 联合概率向量
  • y, 概率向量,必须是随机变量y的概率分布
  • unit,对数化的单位,默认为log2

函数使用:


> x3<- 1:10/sum(1:10)
> y3<- 30:40/sum(30:40)

# 计算条件熵
> CE(x3, y3)
[1] -0.3498852

3.3 JE()联合熵

JE()函数,基于给定的联合概率向量P(X,Y)计算香农的联合熵H(X,Y)。

JE()函数定义:

JE (x, unit = "log2") 

参数列表:

  • x, 联合概率向量
  • unit,对数化的单位,默认为log2

函数使用:

# 创建数据x
> x2 <- 1:100/sum(1:100)

# 联合熵
> JE(x2)
[1] 6.372236

3.4 MI()互信息
MI()函数,根据给定联合概率向量P(X,Y)、概率向量P(X)和概率向量P(X),按公式I(X,Y)= H(X)+ H(Y)-H(X,Y)计算。

函数定义:

MI(x, y, xy, unit = "log2")

参数列表:

  • x, 概率向量
  • x, 概率向量
  • xy, 联合概率向量
  • unit,对数化的单位,默认为log2

函数使用:


# 创建数据集
> x3 <- 1:10/sum(1:10)
> y3<- 20:29/sum(20:29)
> xy3 <- 1:10/sum(1:10)

# 计算互信息
> MI(x3, y3, xy3)
[1] 3.311973

3.5 KL()散度
KL()函数,计算两个概率分布P和Q的Kullback-Leibler散度。
函数定义:

KL(x, test.na = TRUE, unit = "log2", est.prob = NULL)

参数列表:

  • x, 概率向量或数据框
  • test.na, 是否检查NA值
  • unit,对数化的单位,默认为log2
  • est.prob, 用计数向量估计概率的方法,默认值NULL。

函数使用:


# 创建数据集
> df4 <- rbind(x3,y3);df4
         [,1]       [,2]       [,3]       [,4]       [,5]      [,6]      [,7]      [,8]      [,9]
x3 0.01818182 0.03636364 0.05454545 0.07272727 0.09090909 0.1090909 0.1272727 0.1454545 0.1636364
y3 0.08163265 0.08571429 0.08979592 0.09387755 0.09795918 0.1020408 0.1061224 0.1102041 0.1142857
       [,10]
x3 0.1818182
y3 0.1183673

# 计算KL散度 
> KL(df4, unit = "log2") # Default
kullback-leibler 
       0.1392629 
> KL(df4, unit = "log10")
kullback-leibler 
       0.0419223 
> KL(df4, unit = "log")
kullback-leibler 
      0.09652967 

3.5 JSD()散度

JSD()函数,基于具有相等权重的Jensen-Shannon散度,计算距离矩阵或距离值。

公式定义:

函数定义:

JSD(x, test.na = TRUE, unit = "log2", est.prob = NULL)

参数列表:

  • x, 概率向量或数据框
  • test.na, 是否检查NA值
  • unit, 对数化的单位,默认为log2
  • est.prob, 用计数向量估计概率的方法,默认值NULL。

# 创建数据
> x5 <- 1:10
> y5 <- 20:29
> df5 <- rbind(x5,y5)

# 计算JSD
> JSD(df5,unit='log2')
jensen-shannon 
      50.11323 
> JSD(df5,unit='log')
jensen-shannon 
      34.73585 
> JSD(df5,unit='log10')
jensen-shannon 
      15.08559 

# 计算JSD,满足est.prob
> JSD(df5, est.prob = "empirical")
jensen-shannon 
    0.03792749 

3.6 gJSD()散度

gJSD()函数,计算概率矩阵的广义Jensen-Shannon散度。

公式定义:

函数定义:

gJSD(x, unit = "log2", weights = NULL)

参数列表:

  • x, 概率矩阵
  • unit, 对数化的单位,默认为log2
  • weights, 指定x中每个值的权重,默认值NULL。

# 创建数据
> Prob <- rbind(1:10/sum(1:10), 20:29/sum(20:29), 30:39/sum(30:39))

# 计算gJSD
> gJSD(Prob)
[1] 0.023325

4. 应用举例

在我们了解了熵的公式原理和使用方法后,我们就可以做一个案例来试一下。我们定义一个场景的目标:通过用户的看书行为,预测用户是否爱玩游戏。通过我们一步一步地推倒,我们计算出熵,条件熵,联合熵,互信息等指标。

第一步,创建数据集为2列,X列用户看书的类型,包括旅游(Tourism)、美食(Food)、IT技术(IT),Y列用户是否喜欢打游戏,喜欢(Y),不喜欢(N)。


X,Y
Tourism,Y
Food,N
IT,Y
Tourism,N
Tourism,N
IT,Y
Food,N
Tourism,Y

第二步,建立联合概率矩阵,分别计算H(X),Y(X)。

X Y N p(X)
Tourism 2/8=0.25 2/8=0.25 0.25+0.25=0.5
Food 0/8=0 2/8=0.25 0+0.25=0.25
IT 2/8=0.25 0/8=0 0.25+0=0.25
p(Y) 0.25+0+0.25=0.5 0.25+0.25+0=0.5

计算过程


# 分别计算每种情况的概率
p(X=Tourism) = 2/8 + 2/8 = 0.5
p(X=Food) = 2/8 + 0/8 = 0.25
p(X=IT) = 0/8 + 2/8 = 0.25
p(Y=Y) = 4/8 = 0.5
p(Y=N) = 4/8 = 0.5

# 计算H(X)
H(X) = -∑p(xi)*log2(p(xi)) 
 = -p(X=Tourism)*log2(p(X=Tourism) ) -p(X=Food)*log2(p(X=Food) ) -p(X=IT)*log2(p(X=IT) ) 
 = -0.5*log(0.5) -0.25*log(0.25) - 0.25*log(0.25)
 = 1.5

# 计算H(Y)
H(Y) = -∑p(yi)*log2(p(yi)) 
 = -p(Y=Y)*log2(p(Y=Y)) -p(Y=N)*log2(p(Y=N))
 = -0.5*log(0.5) -0.5*log(0.5)
 = 1

第三步,计算每一项的条件熵,H(Y|X=Tourism),H(Y|X=Food),H(Y|X=IT)。


H(Y|X=Tourism) = -p(Y|X=Tourism)*log(p(Y|X=Tourism)) - p(N|X=Tourism)*log(p(N|X=Tourism))
 = -0.5*log(0.5) -0.5*log(0.5)
 = 1

H(Y|X=Food) = -p(Y|X=Food)*log(p(Y|X=Food)) -p(N|X=Food)*log(p(N|X=Food))
 = -0*log(0) -1*log(1)
 = 0

H(Y|X=IT) = -p(Y|X=IT)*log(p(Y|X=IT)) -p(N|X=IT)*log(p(N|X=IT))
 = -1*log(1) -0*log(0) 
 = 0

第四步,计算条件熵H(Y|X)


H(Y|X) = ∑p(xi)*H(Y|xi)
 = p(X=Tourism)*H(Y|X=Tourism) + p(X=Food)*H(Y|X=Food) + p(X=IT)*H(Y|X=IT)
 = 0.5*1 + 0.25*0 + 0.25*0
 = 0.5

第五步,计算联合熵H(X,Y)


H(X,Y) = −∑p(x,y)log(p(x,y))
 = H(X) + H(Y|X)
 = 1.5 + 0.5
 = 2

第六步,计算互信息I(X;Y)


I(X;Y) = H(Y) - H(Y|X)  = 1 - 0.5 = 0.5
 = H(X) + H(Y) - H(X,Y) = 1.5 + 1 - 2 = 0.5

我们把上面的推到过程,用程序来实现一下。


# 创建数据集
> X<-c('Tourism','Food','IT','Tourism','Tourism','IT','Food','Tourism')
> Y<-c('Y','N','Y','N','N','Y','N','Y') 
> df<-cbind(X,Y);df
     X         Y  
[1,] "Tourism" "Y"
[2,] "Food"    "N"
[3,] "IT"      "Y"
[4,] "Tourism" "N"
[5,] "Tourism" "N"
[6,] "IT"      "Y"
[7,] "Food"    "N"
[8,] "Tourism" "Y

变型为频率矩阵


> tf<-table(df[,1],df[,2]);tf
         
          N Y
  Food    2 0
  IT      0 2
  Tourism 2 2

计算概率矩阵


> pX<-margin.table(tf,1)/margin.table(tf);pX
Tourism    Food      IT 
   0.50    0.25    0.25 
> pY<-margin.table(tf,2)/margin.table(tf);pY
  Y   N 
0.5 0.5 
> pXY<-prop.table(tf);pXY
           Y    N
Tourism 0.25 0.25
Food    0.00 0.25
IT      0.25 0.00

计算熵


> H(pX)
[1] 1.5
> H(pY)
[1] 1

# 条件熵 
> CE(pX,pY)
[1] 0.5

# 联合熵 
> JE(pXY)
[1] 2

# 互信息
> MI(pX,pY,pXY)
[1] 0.5

计算原理是复杂的,用R语言的程序实现却是很简单的,几行代码就搞定了,

本文只是对的信息论的初探,重点还是在信息度量方法的R语言实现。信息熵作为信息度量的基本方法,对各种主流的机器学习的算法都有支撑,是我们必须要掌握的知识。了解本质才能发挥数据科学的潜力,学习的路上不断积累和前进。

转载请注明出处:
http://blog.fens.me/r-entropy

打赏作者

R语言实现46种距离算法

R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。

要成为有理想的极客,我们不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让我们一起动起来吧,开始R的极客理想。

关于作者:

  • 张丹(Conan), 程序员/Quant: Java,R,Nodejs
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-distance

前言

距离算法是做数据挖掘常用的一类算法,距离算法有很多种,比如欧式距离、马氏距离、皮尔逊距离,距离算法主要应用在计算数据集之间关系。本文用R语言来philentropy包,实现多种距离的算法,很多可能是大家完全没有听过的,让我们在开拓一下知识领域吧。

目录

  1. 距离算法包philentropy
  2. 46种距离算法详解
  3. 距离函数的使用

1.距离算法包philentropy

在做距离算法调研时,无意中发了philentropy包。它实现了46个不同距离算法和相似性度量,通过不同数据的相似度比较,为基础研究提供了科学基础。philentropy包,为聚类、分类、统计推断、拟合优度、非参数统计、信息理论和机器学习提供了核心的计算框架,支持基于单变量或者多变量的概率函数的计算。

philentropy包主要包括了2种度量的计算方法,距离度量和信息度量。本文介绍距离度量的使用,对于信息度量的使用,请参考文章R语言实现信息度量

philentropy项目github地址:https://github.com/HajkD/philentropy

本文的系统环境为:

  • Win10 64bit
  • R: 3.4.2 x86_64-w64-mingw32

安装philentropy包,非常简单,一条命令就可以了。

~ R
> install.packages("philentropy")
> library(philentropy)

查看距离算法列表

> getDistMethods()
 [1] "euclidean"         "manhattan"         "minkowski"         "chebyshev"        
 [5] "sorensen"          "gower"             "soergel"           "kulczynski_d"     
 [9] "canberra"          "lorentzian"        "intersection"      "non-intersection" 
[13] "wavehedges"        "czekanowski"       "motyka"            "kulczynski_s"     
[17] "tanimoto"          "ruzicka"           "inner_product"     "harmonic_mean"    
[21] "cosine"            "hassebrook"        "jaccard"           "dice"             
[25] "fidelity"          "bhattacharyya"     "hellinger"         "matusita"         
[29] "squared_chord"     "squared_euclidean" "pearson"           "neyman"           
[33] "squared_chi"       "prob_symm"         "divergence"        "clark"            
[37] "additive_symm"     "kullback-leibler"  "jeffreys"          "k_divergence"     
[41] "topsoe"            "jensen-shannon"    "jensen_difference" "taneja"           
[45] "kumar-johnson"     "avg"    

46个距离算法,有一些是我们常用的比如:euclidean,manhattan,minkowski,pearson, cosine,squared_chi, 其他的我也不知道,正好拓宽知识,好好学习一下。

philentropy包的函数,其实很简单,只有14个,大量的算法其实都已经被封装到distance()函数中,直接使用distance()函数就行完成各种算法的计算,让我们使用起来会非常方便。我们来看一下,函数列表:

  • distance(): 计算距离
  • getDistMethods(),获得距离算法列表
  • dist.diversity(),概率密度函数之间的距离差异
  • estimate.probability(),从计数向量估计概率向量
  • lin.cor(),线性相关性判断
  • H(): 香农熵, Shannon’s Entropy H(X)
  • JE() : 联合熵, Joint-Entropy H(X,Y)
  • CE() : 条件熵, Conditional-Entropy H(X|Y)
  • MI() : 互信息, Shannon’s Mutual Information I(X,Y)
  • KL() : KL散度, Kullback–Leibler Divergence
  • JSD() : JS散度,Jensen-Shannon Divergence
  • gJSD() : 通用JS散度,Generalized Jensen-Shannon Divergence
  • binned.kernel.est(),实现了KernSmooth包提供的核密度估计函数的接口

从函数列表来看,主要分为3种类别的函数,第一类是距离测量的函数,包括distance(),
getDistMethods(), dist.diversity(), lin.cor()和 estimate.probability()。第二类是相关性分析,包括lin.cor()函数。第三类是信息度量函数H(),JE(),CE(),MI(),KL(),JSD(),gJSD()。信息度量函数的使用,请参考文章R语言实现信息度量

2. 46种算法详解

接下来,就让我们深入每个算法吧,从名字到公式,再到函数使用,最后到使用场景。

距离算法列表:

  • euclidean:欧式距离,是一个通常采用的距离定义,在m维空间中两个点之间的真实距离,或者向量的自然长度(即该点到原点的距离)。
  • manhattan:曼哈顿距离,用于几何空间度量,表示两个点在标准坐标系上的绝对轴距距离总和。
  • minkowski:闵可夫斯基距离,是欧氏空间中的广义距离函数,其参数p值的不同代表着对空间不同的度量。
  • chebyshev:切比雪夫距离,是向量空间中的一种度量,二个点之间的距离定义是其各坐标数值差绝对值的最大值。
  • sorensen:测量每个样本单位,对单位总数的距离测量贡献度,广告用于生态学。
  • gower:高尔距离,将向量空间缩放为规范化空间,可计算逻辑值,数字,文本的距离,距离结果为0到1之间的数字。
  • soergel:测量每个样本单位,对最大值总数的距离测量贡献度。
  • kulczynski:与soergel相反,测量每个样本单位,对最小值总数的距离测量贡献度。
  • canberra:堪培拉距离,是矢量空间中的点对之间的距离的数值度量,它是L_1距离的加权版本。
  • lorentzian:洛伦兹距离,绝对的差异并应用自然对数。
  • intersection:交叉距离,最小轨道交叉距离,是天文学中用于评估天文物体之间潜在的近距离接近和碰撞风险的度量。它被定义为两个物体的密切轨道的最近点之间的距离。
  • non-intersection:非交叉距离
  • wavehedges:波浪距离,
  • czekanowski:
  • motyka:莫蒂卡方程,是czekanowski的一半。
  • kulczynski_s:
  • tanimoto:是标准化内积的另一种变体。
  • ruzicka:
  • inner_product:内部产品空间,计算两个向量的内积产生标量,有时称为标量积或点积。
  • harmonic_mean:调和平均值。
  • cosine:余弦距离,是用向量空间中两个向量夹角的余弦值作为衡量两个个体间差异的大小的度量。
  • hassebrook(PCE):利用P•Q来测量能量的峰值,简称PCE。
  • jaccard:杰卡德距离,用于计算样本间的相似度,分子是A和B的交集大小,分母是A和B的并集大小。
  • dice:骰子
  • fidelity:保真度,在量子信息理论中,保真度是两个量子态“接近”的度量。
  • bhattacharyya:巴氏距离,测量两个概率分布的相似性。
  • hellinger:海林格,用来度量两个概率分布的相似度,它是F散度的一种。
  • matusita:
  • squared_chord:
  • squared_euclidean:欧式距离的平方
  • pearson:皮尔森距离,分子是两个集合的交集大小,分母是两个集合大小的几何平均值,是余弦距离的一种变型。
  • neyman:奈曼,
  • squared_chi:
  • prob_symm:
  • divergence:散度,
  • clark:克拉克,
  • additive_symm:算术和几何平均散度.
  • kullback-leibler:KL散度,用于计算相熵或信息偏差,是衡量两个分布(P、Q)之间的距离,越小越相似。
  • jeffreys:杰弗里斯,J分歧。
  • k_divergence:K散度,
  • topsoe:托普索,是k_divergence加法的对称形式。
  • jensen-shannon:詹森香农,是topsoe距离的一半。
  • jensen_difference:
  • taneja:塔内加,计算算术和几何平均偏差。
  • kumar-johnson:库马尔-约翰逊,
  • avg:平均

由于精力和基础知识有限,对每一种算法还没有更深入的理解和使用,后面会继续补充。这些距离的详细解释,请参考文章 http://csis.pace.edu/ctappert/dps/d861-12/session4-p2.pdf

这么多种的距离算法,其实可以分成8大距离家族,每个家族中不同的算法思路是类似的,可以通过变形或参数不同赋值,进行算法的相互转换。

L_p Minkowski家族,通过对Minkowski 算法p值的不同赋值,可以转换成不同的算法,当p=1时Minkowski距离转为曼哈顿距离;当p=2变Minkowski距离转为欧氏距离;当p接近极限最大值时,Minkowski距离是转为切比雪夫距离。

L_1家族,用于准确的测量绝对差异的特征。

Intersection 交叉距离家族,用于交叉点之间的相似度变换。

Inter Product 家族,几何空间的相似性度量,用于特定的 P•Q 变量来计算。

Squared-chord 家族,在量子信息理论中,量子态“接近”的度量。

Squared L_2 家族( X^2 Squared 家族),以平方的欧几里得距离做为被除数。

香浓信息熵家族,信息熵偏差测量。

组合公式,利用多种算法思路,进行组合的距离测量方法。

3. 距离函数的使用

了解了这么多的距离算法后,让我们来使用一下philentropy包强大的功能函数,把算法落地。

3.1 distance()函数使用
distance()函数,用来计算两个概率密度函数之间的距离和相似度,上面所列出的所有的距离算法都被封装在了这个函数里。

distance()函数定义:

distance(x, method = "euclidean", p = NULL, test.na = TRUE, unit = "log", est.prob = NULL)

参数列表:

  • x, 数值类型的向量或数据集
  • method, 算法的名称
  • p, minkowski闵可夫斯基距离的p值,p=1为曼哈顿距离,p=2为欧氏距离,p取极限时是切比雪夫距离
  • test.na, 检测数据集是否有NA值,不检测为FALSE,计算会快。
  • unit,对数化的单位,依赖于日志计算的距离
  • est.prob 从计数估计概率,默认值为NULL

计算euclidean距离,用iris的数据集。

> library(magrittr)

# 查看iris数据集
> head(iris)
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
6          5.4         3.9          1.7         0.4  setosa

计算第1个点(第1行)和第2个点(第2行)的euclidean距离,分别使用philentropy包的distance(),Stats包的dist(),和自己通过公式计算。

# 使用distance()函数
> dat1<-iris[1:2,-5]
> distance(dat1, method="euclidean")
Metric: 'euclidean' using unit: 'log'.
euclidean 
0.5385165 

# 再使用系统自带的dist()函数
> dist(dat1)
          1
2 0.5385165

# 公式计算
> (dat1[1,]-dat1[2,])^2 %>% sum %>% sqrt
[1] 0.5385165

3种方法,计算的结果是完全一致。

接下来,我们构建一个iris的距离矩阵,为了展示清楚,我们选iris的前6个点来计算距离。分别使用distance()和dist()函数。


> dat2<-head(iris[,-5])

# 距离矩阵
> distance(dat2)
Metric: 'euclidean' using unit: 'log'.
          v1        v2       v3        v4        v5        v6
v1 0.0000000 0.5385165 0.509902 0.6480741 0.1414214 0.6164414
v2 0.5385165 0.0000000 0.300000 0.3316625 0.6082763 1.0908712
v3 0.5099020 0.3000000 0.000000 0.2449490 0.5099020 1.0862780
v4 0.6480741 0.3316625 0.244949 0.0000000 0.6480741 1.1661904
v5 0.1414214 0.6082763 0.509902 0.6480741 0.0000000 0.6164414
v6 0.6164414 1.0908712 1.086278 1.1661904 0.6164414 0.0000000

# 下三角距离矩阵
> dist(dat2)
          1         2         3         4         5
2 0.5385165                                        
3 0.5099020 0.3000000                              
4 0.6480741 0.3316625 0.2449490                    
5 0.1414214 0.6082763 0.5099020 0.6480741          
6 0.6164414 1.0908712 1.0862780 1.1661904 0.6164414

验证后,我们就可以放心使用distance()函数。通过对比实验,我们可以很快的学习并使用各种距离算法。

3.2 dist.diversity()函数
dist.diversity()函数,用来计算所有距离的值。由于有一些距离有对于数据集本身的要求,所以我们需要构建一个能适应所有距离算法的数据集。


# 生成数据集,2个点,10个维度
> P <- 1:10/sum(1:10)
> Q <- 20:29/sum(20:29)
> x <- rbind(P,Q)

# 打印数据集
> head(x)
        [,1]       [,2]       [,3]       [,4]       [,5]      [,6]      [,7]      [,8]
P 0.01818182 0.03636364 0.05454545 0.07272727 0.09090909 0.1090909 0.1272727 0.1454545
Q 0.08163265 0.08571429 0.08979592 0.09387755 0.09795918 0.1020408 0.1061224 0.1102041
       [,9]     [,10]
P 0.1636364 0.1818182
Q 0.1142857 0.1183673

使用dist.diversity()函数计算所有的距离。

> dist.diversity(x,p=2)
euclidean         manhattan         minkowski         chebyshev          sorensen 
       0.12807130        0.35250464        0.12807130        0.06345083        0.17625232 
            gower           soergel      kulczynski_d          canberra        lorentzian 
       0.03525046        0.29968454        0.42792793        2.09927095        0.34457827 
     intersection  non-intersection        wavehedges       czekanowski            motyka 
       0.82374768        0.17625232        3.16657887        0.17625232        0.58812616 
     kulczynski_s          tanimoto           ruzicka     inner_product     harmonic_mean 
       2.33684211        0.29968454        0.70031546        0.10612245        0.94948528 
           cosine        hassebrook           jaccard              dice          fidelity 
       0.93427641        0.86613103        0.13386897        0.07173611        0.97312397 
    bhattacharyya         hellinger          matusita     squared_chord squared_euclidean 
       0.02724379        0.32787819        0.23184489        0.05375205        0.01640226 
          pearson            neyman       squared_chi         prob_symm        divergence 
       0.16814418        0.36742465        0.10102943        0.20205886        1.49843905 
            clark     additive_symm  kullback-leibler          jeffreys      k_divergence 
       0.86557468        0.53556883        0.09652967        0.22015096        0.02922498 
           topsoe    jensen-shannon jensen_difference            taneja     kumar-johnson 
       0.05257867        0.02628933        0.02628933        0.02874841        0.62779644 
              avg 
       0.20797774 

3.3 estimate.probability()函数
estimate.probability()函数,采用数字计数来计算向量的估计概率。estimate.probability()函数方法实现,目前只有一个方法实现,计算每个值占合计的比率,其实就是一种数据标准归的计算方法。

> estimate.probability
function (x, method = "empirical") 
{
    if (!is.element(method, c("empirical"))) 
        stop("Please choose a valid probability estimation method.")
    if (method == "empirical") {
        return(x/sum(x))
    }
}
>environment: namespace:philentropy<

我们新建一个向量,用estimate.probability()函数,来计算向量的估计概率。

# 新建x1向量
> x1<-runif(100);head(x1)
[1] 0.6598775 0.2588441 0.5329965 0.5294842 0.8331355 0.3326702

# 计算估计概率
> x2<-estimate.probability(x1);head(x2)
[1] 0.013828675 0.005424447 0.011169702 0.011096097 0.017459543 0.006971580

# 打印统计概率
> summary(x2)
Min.   1st Qu.    Median      Mean   3rd Qu.      Max.
0.0002181 0.0044882 0.0096318 0.0100000 0.0153569 0.0206782

# 画散点图
> plot(x1,x2)

从图中看到,整个数据分布在对角线上,x1为均匀分布生成的向量,x2为x1的估计概率,仅仅做了做数据值域进行了缩放,并没有影响数据的分布变化和线性特征,其实就是对数据做了一个标准化的过程。

3.4 lin.cor()函数
lin.cor()函数,用来计算两个向量之间的线性相关,或计算矩阵的相关矩阵。

函数定义:

lin.cor(x, y = NULL, method = "pearson", test.na = FALSE)

参数列表:

  • x,变量1
  • y,变量2,与x进行比较
  • method,相关性算法的名称,默认为pearson距离算法,支持5种算法分为是pearson,pearson2,sq_pearson,kendall,spearman
  • test.na, 检测数据集是否有NA值,不检测为FALSE,计算会快。

相关性计算,最大值1为完全正相关,最小值-1为完全负相关,0为不相关。我们来创建数据集,进行相关性的测试。

# 创建向量x1,x2,x3
> x1<-runif(100)
> x2<-estimate.probability(x1)
> x3<-rnorm(100)

# 判断x1,x2的相关性,pearson 皮尔森相关系数
> lin.cor(x1,x2)
pearson 
      1 

# 判断x1,x3的相关性,pearson 皮尔森相关系数
> lin.cor(x1,x3)
  pearson 
0.0852527 

# 判断x1,x3的相关性,pearson2 皮尔森非集中的相关系数
> lin.cor(x1,x3,method = 'pearson2')
  pearson2 
0.01537887 

# 判断x1,x3的相关性,sq_pearson 皮尔森平方的相关系数
> lin.cor(x1,x3,method = 'sq_pearson')
sq_pearson 
0.00151915 

# 判断x1,x3的相关性,kendall 肯德尔相关系数
> lin.cor(x1,x3,method = 'kendall')
kendall 
      0 

# 判断x1,x3的相关性,spearman 斯皮尔曼相关系数
> lin.cor(x1,x3,method = 'spearman')
spearman 
       0 

通过lin.cor()函数,可以快速进行线性相关性的验证,非常方便。

本文重点介绍了philentropy包,对于距离算法的定义和距离测量的函数的使用。很多的距离算法我也是第一次学习,知识需要积累和总结,本文不完善的内容,后面我们找时间再进行补充。如果本文描述有不当的地方,也请各位朋友,给予指点,让我们一起把知识进行积累。

下一篇文章我们将介绍信息度量函数的使用和算法,请大家继续阅读文章用R语言实现信息度量

转载请注明出处:
http://blog.fens.me/r-distance

打赏作者

R语言轻巧的时间包hms

R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。

要成为有理想的极客,我们不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让我们一起动起来吧,开始R的极客理想。

关于作者:

  • 张丹(Conan), 程序员/Quant: Java,R,Nodejs
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-hms/

前言
时间是数据的基本维度,是在做数据处理的时候,必须要掌握的技术。根据时间周期的不同,通常把时间分为,年、月、日、时、分、秒、毫秒等。对于年月日的数据是最常见的,也有很多的处理工具,时分秒的数据通常也会用处理日期的工具,这样有时候就不太方便。

hms包,很小很轻,专注于时、分、秒的时间数据处理。

目录

  1. hms包介绍
  2. hms包的使用

1. hms包介绍

hms包,用于存储和格式化时间,基于difftime类型,使用S3的面向对象数据结构。

本文的系统环境为:

  • Win10 64bit
  • R: 3.4.2 x86_64-w64-mingw32

安装hms包,非常简单,一条命令就可以了。


~ R
> install.packages("hms")
> library(hms)

函数列表:

  • hms: 创建一个hms类型对象
  • is.hms: 判断是否是hms类型
  • parse_hm: 解析hm值
  • parse_hms: 解析hms值
  • round_hms:四舍五入对齐
  • trunc_hms:裁剪对齐
  • as.hms: hms泛型函数,S3类型,用于隐式调用它的继承函数
  • as.hms.character: character转型hms,用于as.hms的继承调用
  • as.hms.default: hms转型hms,用于as.hms的继承调用
  • as.hms.difftime:difftime转型hms,用于as.hms的继承调用
  • as.hms.numeric: numeric转型hms,用于as.hms的继承调用
  • as.hms.POSIXlt: POSIXlt转型hms,用于as.hms的继承调用
  • as.hms.POSIXt: POSIXt转型hms,用于as.hms的继承调用
  • as.character.hms: hms转型character,用于as.character的继承调用
  • as.data.frame.hms: hms转型data.frame,用于as.data.frame的继承调用
  • as.POSIXct.hms: hms转型POSIXct,用于as.POSIXct的继承调用
  • as.POSIXlt.hms: hms转型POSIXlt,用于as.POSIXlt的继承调用
  • format.hms: 格式化hms,用于format的继承调用
  • print.hms: 打印hms对像,用于print的继承调用

从函数列表可以看到,hms包的功能很单一,就在做数据类型和数据变型,是底层的数据结构包,设计思路与zoo包的设计思路一致。zoo包的详细介绍,请参考文章R语言时间序列基础库zoo

hms包中,有大量的as.xxx()函数、format.hms()函数和print.hms()函数,是被用于S3类型函数继承调用的,是不需要我们在写程序的时候显示调用的。S3数据结构详细介绍,请参考文章R语言基于S3的面向对象编程

2. hms包的使用

接下来,我们找几个重点函数进行介绍。

2.1 hms()函数
hms()函数,用于创建一个hms类型的对象。
函数定义:

hms(seconds = NULL, minutes = NULL, hours = NULL, days = NULL)

hms()函数,接收4个参数,分别对应秒,分,时,日。

创建hms对象


# 创建12:34:56的时间对象
> a1<-hms(56, 34, 12);a1
12:34:56

# 创建 10日12:34:56的时间对象
> a2<-hms(56, 34, 12,10);a2
252:34:56

打印结果的第一位252=10*24+12。

2.2 is.hms: 判断是否是hms类型


# 判断是否是hms类型
> is.hms(a1)
[1] TRUE

# 查看hms类型,父类是difftime类型
> class(a1)
[1] "hms"      "difftime"

# 查看hms的属性
> attributes(a1)
$units
[1] "secs"

$class
[1] "hms"      "difftime"

# 查看hms对象的静态数据结构
> str(a1)
Classes 'hms', 'difftime'  atomic [1:1] 45296
  ..- attr(*, "units")= chr "secs"

# 查看面向对象类型
> library(pryr)
> otype(a1)
[1] "S3"

2.3 as.xxx.hms:把hms转型到其他类型


# 默认转型
> as.hms(a1)
12:34:56

# hms转型character,实际会隐式调用as.character.hms()函数
> as.character(a1)
[1] "12:34:56"

# hms转型POSIXct
> as.POSIXct(a1)
[1] "1970-01-01 12:34:56 UTC"

# hms转型POSIXlt
> as.POSIXlt(a1)
[1] "1970-01-01 12:34:56 UTC"

由于我们没有定义as.Date.hms()函数,所以as.Date()函数,不能认识hms类型的转换。


# hms转型Date
> as.Date(a1)
Error in as.Date.default(a1) : 不知如何将'a1'转换成“Date”类别
Error during wrapup: cannot open the connection

自己定义一个as.Date.hms()函数,仅用于转型实验,没有实际业务意义。


# 函数定义
> as.Date.hms<-function(hms){
+   s<-paste(Sys.Date(),' ',hms,sep="")
+   as.Date(s)
+ }

# 显示调用函数
> as.Date.hms(a1)
[1] "2018-12-14"

# 隐式调用函数
> as.Date(a1)
[1] "2018-12-14"

2.4 as.hms.xxx:把其他类型转型到hms


# 把字符串转hms
> as.hms('19:13:14')
19:13:14
# 非法时间字符串转型
> as.hms('19:78:14')
NA

# 数字转型
> as.hms(111312)
30:55:12

# 时间转型
> as.hms(Sys.time())
14:22:59.462795

# 日期转型,同样发生了错误
> as.hms(Sys.Date())
Error: Can't convert object of class Date to hms.
Error during wrapup: cannot open the connection

2.5 parse_hms()/parse_hm()字符串解析

parse_hms对字符串进行转型,对比parse_hms()与as.hms()结果一样的。


# 执行parse_hms
> parse_hms("12:34:56.789")
12:34:56.789
> as.hms("12:34:56.789")
12:34:56.789

# 执行parse_hm
> parse_hm("12:34")
12:34:00
> as.hms("12:34")
NA

打印parse_hms 函数名,查看源代码实现。


> parse_hms 
function (x) {
as.hms(as.difftime(as.character(x), format = "%H:%M:%OS",
units = "secs"))
}
>environment: namespace:hms<

parse_hms()函数,实际就是调用了as.hms()函数。

2.6 round_hms/trunc_hms
round_hms()函数,是把时间进行四舍五入对齐。


# 按秒,以5的倍数进行对齐,四舍五入
> round_hms(as.hms("12:34:51"), 5)
12:34:50
> round_hms(as.hms("12:34:54"), 5)
12:34:55
> round_hms(as.hms("12:34:56"), 5)
12:34:55
> round_hms(as.hms("12:34:59"), 5)
12:35:00

# 按秒,以60的倍数对齐
> round_hms(as.hms("12:34:56"), 60)
12:35:00

trunc_hms()函数,是把时间进行裁剪对齐。


# 按秒去掉末位,以5的倍数进行对齐
> trunc_hms(as.hms("12:34:01"), 5)
12:34:00
> trunc_hms(as.hms("12:34:44"), 5)
12:34:40
> trunc_hms(as.hms("12:34:56"), 60)
12:34:00

2.7 在data.frame中插入hms列


# 创建data.frame
> df<-data.frame(hours = 1:3, hms = hms(hours = 1:3))
> df
  hours      hms
1     1 01:00:00
2     2 02:00:00
3     3 03:00:00

# 查看df的静态结构
> str(df)
'data.frame':	3 obs. of  2 variables:
 $ hours: int  1 2 3
 $ hms  :Classes 'hms', 'difftime'  atomic [1:3] 3600 7200 10800
  .. ..- attr(*, "units")= chr "secs"

hms包很轻巧很简单,但却可以快速提高帮助我们处理时分秒数据,这些基础函数库是需要我们完全掌握和熟练运用的。

转载请注明出处:
http://blog.fens.me/r-hms/

打赏作者

用R语言实现密度聚类dbscan

R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。

要成为有理想的极客,我们不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让我们一起动起来吧,开始R的极客理想。

关于作者:

  • 张丹(Conan), 程序员/Quant: Java,R,Nodejs
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-cluster-dbscan

前言

聚类是一种将数据点按一定规则分群的机器学习技术,k-Means聚类是被用的最广泛的也最容易理解的一种。除了K-Means的方法,其实还有很多种不同的聚类方法,本文将给大家介绍基于密度的聚类,我们可以通过使用dbscan包来实现。

目录

  1. DBSCAN基于密度的聚类
  2. dbscan包介绍
  3. kNN()函数使用
  4. dbscan()函数使用
  5. hdbscan()函数使用

1. DBSCAN基于密度的聚类

DBSCAN(Density-Based Spatial Clustering of Applications with Noise)聚类算法,它是一种基于高密度连通区域的基于密度的聚类算法,能够将具有足够高密度的区域划分为簇,并在具有噪声的数据中发现任意形状的簇。

DBSCAN需要两个重要参数:epsilon(eps)和最小点(minPts)。参数eps定义了点x附近的邻域半径ε,它被称为x的最邻居。参数minPts是eps半径内的最小邻居数。

上图中(a),数据集中的任何点x邻居(6=minPts)都被标记为核心点,ε是半径。上图中(b),x为核心点,y的邻居小于(4<minpts)是边界点,但它属于核心点x的最邻居。z点既不是核心也不是边界点,它被称为噪声点或异常值。

dbscan算法将数据点分为三类:

  • 核心点:在半径eps内含有超过minPts数目的点。
  • 边界点:在半径eps内点的数量小于使用DBSCAN进行聚类的时候,不需要预先指定簇的个数,最终的簇的个数不确定。minPts,但是落在核心点的邻域内的点。
  • 噪音点:既不是核心点也不是边界点的点

DBSCAN算法的执行过程

1、DBSCAN算法随机从一个未被访问的数据点x开始,以eps为半径搜索范围内的所有邻域点。

2、如果x点在该邻域内有足够数量的点,数量大于等于minPts,则聚类过程开始,并且当前数据点成为新簇中的第一个核心点。否则,该点将被标记为噪声。该点都会被标记为“已访问”。

3、新簇中的每个核心点x,它的eps距离邻域内的点会归为同簇。eps邻域内的所有点都属于同一个簇,然后对才添加到簇中的所有新点重复上述过程。

4、重复步骤2和3两个过程,直到确定了簇中的所有点才停止,即访问和标记了聚类的eps邻域内的所有点。

5、当完成了这个簇的划分,就开始处理新的未访问的点,发现新的簇或者是噪声。重复上述过程,直到所有点被标记为已访问才停止。这样就完成了,对所有点的聚类过程。

优点和缺点

DBSCAN具有很多优点,提前不需要确定簇的数量。不同于Mean-shift算法,当数据点非常不同时,会将它们单纯地引入簇中,DBSCAN能将异常值识别为噪声。另外,它能够很好地找到任意大小和任意形状的簇。

DBSCAN算法的主要缺点是,当数据簇密度不均匀时,它的效果不如其他算法好。这是因为当密度变化时,用于识别邻近点的距离阈值ε和minPoints的设置将随着簇而变化。在处理高维数据时也会出现这种缺点,因为难以估计距离阈值eps。

2. dbscan包介绍

dbscan包,提供了基于密度的有噪声聚类算法的快速实现,包括 DBSCAN(基于密度的具有噪声的应用的空间聚类),OPTICS(用于识别聚类结构的排序点),HDBSCAN(分层DBSCAN)和LOF(局部异常因子)算法,dbscan底层使用C++编程,并建立kd树的数据结构进行更快的K最近邻搜索,从而实现加速。

本文的系统环境为:

  • Win10 64bit
  • R 3.4.2 x86_64

dbscan包的安装非常简单,只需要一条命令就能完成。


~ R
> install.packages("dbscan")
> library(dbscan)

函数列表:

  • dbscan(), 实现DBSCAN算法
  • optics(), 实现OPTICS算法
  • hdbscan(), 实现带层次DBSCAN算法
  • sNNclust(), 实现共享聚类算法
  • jpclust(), Jarvis-Patrick聚类算法
  • lof(), 局部异常因子得分算法
  • extractFOSC(),集群优选框架,可以通过参数化来执行聚类。
  • frNN(), 找到固定半径最近的邻居
  • kNN(), 最近邻算法,找到最近的k个邻居
  • sNN(), 找到最近的共享邻居数量
  • pointdensity(), 计算每个数据点的局部密度
  • kNNdist(),计算最近的k个邻居的距离
  • kNNdistplot(),画图,最近距离
  • hullplot(), 画图,集群的凸壳

dbscan包,提供了多个好用的函数,我们接下来先介绍3个函数,分别是kNN(),dbscan(), hdbscan(),其他的函数等以后有时间,再单独进行使用介绍。

3. kNN()函数使用

kNN()函数,使用kd-tree数据结构,用来快速查找数据集中的所有k个最近邻居。

函数定义:


kNN(x, k, sort = TRUE, search = "kdtree", bucketSize = 10, splitRule = "suggest", approx = 0)

参数列表

  • x,数据矩阵,dist对象或kNN对象。
  • k,要查找的邻居数量。
  • sort,按距离对邻居进行排序。
  • search,最近邻搜索策略,使用kdtree,linear或dist三选一,默认为kdtree。
  • bucketSize,kd-tree叶子节点的最大值。
  • splitRule,kd-tree的拆分规则,默认用SUGGEST。
  • approx,使用近似方法,加速计算。

函数使用:以iris鸢尾花的数据集,做为样本。聚类是不需要有事前有定义的,所以我们把iris的种属列去掉。

# 去掉种属列
> iris2 <- iris[, -5]
> head(iris2)
  Sepal.Length Sepal.Width Petal.Length Petal.Width
1          5.1         3.5          1.4         0.2
2          4.9         3.0          1.4         0.2
3          4.7         3.2          1.3         0.2
4          4.6         3.1          1.5         0.2
5          5.0         3.6          1.4         0.2
6          5.4         3.9          1.7         0.4

使用kNN()函数,来计算iris2数据集中,每个值最近的5个点。


# 查询最近邻的5个点
> nn <- kNN(iris2, k=5)

# 打印nn对象
> nn
k-nearest neighbors for 150 objects (k=5).
Available fields: dist, id, k, sort

# 查询nn的属性列表
> attributes(nn)
$names
[1] "dist" "id"   "k"    "sort"

$class
[1] "kNN" "NN" 

打印出,每个点最近邻的5个点。行,为每个点索引值,列,为最近邻的5个点,输出的矩阵为索引值。


> head(nn$id)
      1  2  3  4  5
[1,] 18  5 40 28 29
[2,] 35 46 13 10 26
[3,] 48  4  7 13 46
[4,] 48 30 31  3 46
[5,] 38  1 18 41  8
[6,] 19 11 49 45 20

打印出,每个点与最近的5个点的距离值。行,为每个点的索引,列,为最近邻的5个点,输出的矩阵为距离值。


> head(nn$dist)
             1         2         3         4         5
[1,] 0.1000000 0.1414214 0.1414214 0.1414214 0.1414214
[2,] 0.1414214 0.1414214 0.1414214 0.1732051 0.2236068
[3,] 0.1414214 0.2449490 0.2645751 0.2645751 0.2645751
[4,] 0.1414214 0.1732051 0.2236068 0.2449490 0.2645751
[5,] 0.1414214 0.1414214 0.1732051 0.1732051 0.2236068
[6,] 0.3316625 0.3464102 0.3605551 0.3741657 0.3872983

如果我们要查看索引为33的点,与哪5个点最紧邻,可以用下面的方法。


# 设置索引
> idx<-33

# 打印与33,最近邻的5个点的索引
> nn$id[idx,]
 1  2  3  4  5 
34 47 20 49 11 

# 画图
> cols = ifelse(1:nrow(iris2) %in% nn$id[idx,],"red", "black")
> cols[idx]<-'blue'
> plot(iris2,pch = 19, col = cols)

我们的数据集是多列的,把每2列组合形成的二维平面,都进行输出。蓝色表示索引为33的点,红色表示最紧邻的5个点,黑色表示其他的点。

从图中,可以很直观的看到,这几点确实是密集的在一起,也就是找到了最近邻。

接下来,我们画出连线图,选取第一列(Sepal.Length)和第二列(Sepal.Width),按取画出最紧邻前5连接路径。

> plot(nn, iris2)

通过连接路径,我们就能很清晰的看到,最紧邻算法的分组过程,连接在一起的就够成了一个分组,没有连接在一起的就是另外的分组,上图中可以看出来分成了2个组。

再对nn进行二次最近邻计算,画出前2的连接路径。

> plot(kNN(nn, k = 2), iris2)

通过2次的最紧邻缩减,连接路径大幅度减少了,又形成了新的独立区块。

2. dbscan()函数使用

dbscan是一种基于密度的聚类算法,这类密度聚类算法一般假定类别可以通过样本分布的紧密程度决定。同一类别的样本,他们之间的紧密相连的,也就是说,在该类别任意样本周围不远处一定有同类别的样本存在。

函数定义:

dbscan(x, eps, minPts = 5, weights = NULL, borderPoints = TRUE, ...)

参数解释:

  • x, 矩阵或者距离对象,frNN对象。
  • eps,半径的大小。
  • minPts, 半径区域中的最小点数量,默认为5
  • weights, 数据点的权重,仅用于加权聚类
  • borderPoints,边界点是否为噪声,默认为TRUE;为FALSE时,边界点为噪声。
  • …,将附加参数传递给固定半径最近邻搜索算法,调用frNN。

函数使用:以iris鸢尾花的数据集,做为样本。聚类是不需要有事前有定义的,所以我们把iris的种属列去掉。

# 去掉种属列
> iris2 <- iris[, -5]
> head(iris2)
  Sepal.Length Sepal.Width Petal.Length Petal.Width
1          5.1         3.5          1.4         0.2
2          4.9         3.0          1.4         0.2
3          4.7         3.2          1.3         0.2
4          4.6         3.1          1.5         0.2
5          5.0         3.6          1.4         0.2
6          5.4         3.9          1.7         0.4

在使用dbscan函数时,我们要输出2个参数,eps和minPts。

  • eps,值可以使用绘制k-距离曲线(k-distance graph)方法得到,在k-距离曲线图明显拐点位置为较好的参数。若参数设置过小,大部分数据不能聚类;若参数设置过大,多个簇和大部分对象会归并到同一个簇中。
  • minPts,通常让minPts≥dim+1,其中dim表示数据集聚类数据的维度。若该值选取过小,则稀疏簇中结果由于密度小于minPts,从而被认为是边界点儿不被用于在类的进一步扩展;若该值过大,则密度较大的两个邻近簇可能被合并为同一簇。

下面我们通过绘制k-距离曲线,寻找knee,即明显拐点位置为对应较好的参数,找到适合的eps值。使用kNNdistplot()函数,让参数k=dim + 1,dim为数据集列的个数,iris2是4列,那么设置k=5。

# 画出最近距离图
> kNNdistplot(iris2, k = 5)
> abline(h=0.5, col = "red", lty=2)

kNNdistplot()会计算点矩阵中的k=5的最近邻的距离,然后按距离从小到大排序后,以图形进行展示。x轴为距离的序号,y轴为距离的值。图中黑色的线,从左到右y值越来越大。

通过人眼识别,k-距离曲线上有明显拐点,我们以y=0.5平行于x轴画一条红色线,突出标识。所以,最后确认的eps为0.5。

调用dbscan()函数,进行对iris2数据集进行聚类,eps=0.5,minPts=5。

> res <- dbscan(iris2, eps = 0.5, minPts = 5)
> res
DBSCAN clustering for 150 objects.
Parameters: eps = 0.5, minPts = 5
The clustering contains 2 cluster(s) and 17 noise points.

 0  1  2 
17 49 84 

Available fields: cluster, eps, minPts

聚类后,一共分成了2组,第1组49个值,第2组84个值,另外,第0组17个值为噪声点。把聚类的结果画图展示。

> pairs(iris, col = res$cluster + 1L)

数据集是多列的,把每2列组合形成的二维平面,都进行输出。红色点表示第1组,绿色点表示为第2组,黑色点表示噪声点。这样就完成了有噪声的基于密度的dbscan聚类。

5. hdbscan()函数使用

hdbscan(),快速实现了分层DBSCAN算法,与stats包中的hclust()方法形成的传统分层聚类方法类似。

函数定义:

hdbscan(x, minPts, xdist = NULL,gen_hdbscan_tree = FALSE, gen_simplified_tree = FALSE)

参数解释:

  • x,矩阵或者距离对象
  • minPts,区域中的最小点数量
  • xdist,dist对象,可以提前算出来,当参数传入
  • gen_hdbscan_tree,生成一个hdbscan树
  • gen_simplified_tree,生成一个简化的树结构

5.1 iris鸢尾花的数据集
以iris鸢尾花的数据集,做为样本,去掉种属列。设置minPts =5让当前群集中最小的数量为5,开始聚类。

> hcl<-hdbscan(iris2, minPts = 5);hcl
HDBSCAN clustering for 150 objects.
Parameters: minPts = 5
The clustering contains 2 cluster(s) and 0 noise points.

  1   2 
100  50 

Available fields: cluster, minPts, cluster_scores, membership_prob, outlier_scores, hc

聚类后,一共分成了2组,第1组100个值,第2组50个值,没有噪声点。生成的hcl对象包括6个属性。
属性解释

  • cluster,表明属性哪个群集,零表示噪声点。
  • minPts,群集中最小的数量
  • cluster_scores,每个突出(“平坦”)群集的稳定性分数之和。
  • membership_prob,群集内某点的“概率”或个体稳定性
  • outlier_scores,每个点的异常值
  • hc,层次结构对象

把聚类的结果画图展示。

> plot(iris2, col=hcl$cluster+1, pch=20)

数据集是多列的,把每2列组合形成的二维平面,都进行输出。红色点表示第1组,绿色点表示为第2组,这样就完成了hdbscan聚类。

打印hcl对象层次结构,包括150个数据,聚法方法是健壮单一的,距离是相互可达。

> hcl$hc

Call:
hdbscan(x = iris2, minPts = 5)

Cluster method   : robust single 
Distance         : mutual reachability 
Number of objects: 150 

画出层次的合并过程图

> plot(hcl$hc, main="HDBSCAN* Hierarchy")

从图可以清楚的看出,主要的2类的分支,区分度比较高。

5.2 moons数据集
由于iris数据集用hdbscan聚类获得的结果,与真实的数据分类结果不一致。我们再用dbscan包自带的数据集moons做一下测试。

先准备数据,加载moons数据集,了解数据基本情况,画出散点图。

# 加载dbscan自带数据集
> data("moons")
> head(moons)
            X          Y
1 -0.41520756  1.0357347
2  0.05878098  0.3043343
3  1.10937860 -0.5097378
4  1.54094828 -0.4275496
5  0.92909498 -0.5323878
6 -0.86932470  0.5471548

# 画出散点图
> plot(moons, pch=20)

用hdbscan()函数,实现层次dbscan算法。


> cl <- hdbscan(moons, minPts = 5)
> cl
HDBSCAN clustering for 100 objects.
Parameters: minPts = 5
The clustering contains 3 cluster(s) and 0 noise points.

 1  2  3 
25 25 50 

Available fields: cluster, minPts, cluster_scores, membership_prob, outlier_scores, hc

一共100条数据,被分成了3类,没有噪声。把聚类的结果画图展示。


# 画图
> plot(moons, col=cl$cluster+1, pch=20)

打印层次结构


> cl$hc
Call:
hdbscan(x = moons, minPts = 5)

Cluster method   : robust single 
Distance         : mutual reachability 
Number of objects: 100 

画出层次的合并过程图

> plot(cl$hc, main="HDBSCAN* Hierarchy")

从图可以清楚的看出,主要的3类的分支,区分度比较高。

如果我们想省略分层的细节,我们可以只画出主要分支,并标识类别。

plot(cl, gradient = c("purple", "blue", "green", "yellow"), show_flat = T)

接下来,我们要对群集的稳定性做一些优化,cluster_scores属性可以查看集群的得分。

> cl$cluster_scores
        1         2         3 
110.70613  90.86559  45.62762 

通过membership_prob属性,画图表示个体的稳定性。

# 打印membership_prob
> head(cl$membership_prob)
[1] 0.4354753 0.2893287 0.4778663 0.4035933 0.4574012 0.4904582

# 计算群集的数量
> num<-length(cl$cluster_scores)

# 从彩虹色中取得对应数量的颜色
> rains<-rainbow(num)
> cols<-cl$cluster
> cols[which(cols==1)]<-rains[1]
> cols[which(cols==2)]<-rains[2]
> cols[which(cols==3)]<-rains[3]

# 设置透明度,表示个体的稳定性
> plot(moons, col=alpha(cols,cl$membership_prob), pch=19)

最后,我们可以在图中,在标记出异常值得分最高的前6个点。

# 对异常值进行排序,取得分最高的
> top_outliers <- order(cl$outlier_scores, decreasing = TRUE) %>% head
> plot(moons, col=alpha(cols,cl$outlier_scores), pch=19)
> text(moons[top_outliers, ], labels = top_outliers, pos=3)

从图中看到,异常得分高的点(outlier_scores)与个体的稳定性(membership_prob),并不是同一类点。异常值通常被认为是,偏离其假定的基础分布的离群点。

通过上面3个函数的使用案例,我们了解了如何用dbscan包实现基于密度的聚类方法。真实世界的数据是复杂的,我们用来分析数据的工具也是多样的,多掌握一种工具、多一些知识积累,让我们迎接真实世界数据的挑战吧。

转载请注明出处:
http://blog.fens.me/r-cluster-dbscan

打赏作者