• Archive by category "R语言实践"
  • (Page 2)

Blog Archives

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

打赏作者

2018南京信息工程大学:R语言与计量经济学

跨界知识聚会系列文章,“知识是用来分享和传承的”,各种会议、论坛、沙龙都是分享知识的绝佳场所。我也有幸作为演讲嘉宾参加了一些国内的大型会议,向大家展示我所做的一些成果。从听众到演讲感觉是不一样的,把知识分享出来,你才能收获更多。

关于作者

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

转载请注明出处:
http://blog.fens.me/meeting-nuist-20181201

前言

学生面临毕业的竞争,如何能脱颖而了,不仅需要有良好的知识储备,更要具有动手能力,同时了解市场的需求是什么。很多同学,在学校读书学习,最后都变成了,“死读书,读死书”,完成不知道社会是什么样的,只是把书本上和老师教的内容掌握,这是完全不行的!!!

当老师充分的认识到,学生需要更多的知识来源时,就会选择与行业专家来合作,从而打开学生的视野。那么,本次的教学就是这样的一个目的,南京信息工程大学经济管理学院的鲁老师找到我,希望我给他的学生带来一次,R语言与计量经济的培训,让学生不仅学到知识,更多看看外面的世界是什么样子的。

目录

  1. 我的培训主题:基于R语言的数据挖掘与计量经济学
  2. 会议体验和照片分享

1. 我的培训主题:基于R语言的数据挖掘与计量经济学

50人的小班课程,定位就是把R语言和计量经济学的知识结合起来,让学生们提高动手能力,掌握一门技能,提高毕业时求职的竞争力。本次培训的活动主页

课程大纲:2018年全国高校计量经济学及R语言应用师资培训大纲

本期培训内容以“计量经济学及应用”考核大纲为基础,重点讲解R语言编程、金融数据计算、数据分析与可视化等内容,并对2018年度计量经济学及应用试题进行实操演练。

经过2天的密集的授课,让同学们从3个维度获得了提升:

  • IT思维:熟悉的R语言编程
  • 统计思维:把计量经济学的模型通过R语言落地
  • 金融市场思维:了解金融市场,把书里上的知识与市场进行了对接。

2. 会议体验和照片分享

我的介绍和照片分享。

张丹,青萌数海CTO,微软MVP,资深R语言技术专家,,《R的极客理想-量化投资篇》作者,微软MVP。10年编程经验,获得10项SUN及IBM技术认证。前民生银行大数据分析师。个人博客 http://fens.me, Alexa全球排名70K。

2.2 会议体验证和相关照片

本次的场地南京信息工程大学经济管理学院,由学院的鲁老师负责制定整体的教学计划。整体来说,课程安排的很紧很充实,也很累。嗓子疼!!又赶上了南京的大雾霾天,和北京真是有的一比了。

大合照

我在讲课

上课啦

动手实践,讨论问题

教师里问问题的学生

学校正门

体育场

图书馆

据说全国排名1-2的气象学院

野鸭湖

最后,课程获得了,同学们和老师的好评!希望同学们能消化吸收,未来的世界是属于你们的。

转载请注明出处:
http://blog.fens.me/meeting-nuist-20181201

打赏作者

用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-discretization

前言

在做数据挖掘模型的时候,我们有时会需要把连续型变量转型离散变量,这种转换的过程就是数据离散化,分箱就是离散化常用的一种方法。

数据离散化处理属于数据预处理的一个过程,R语言在数据处理上有天然的优势,也有直接用于离散化计算的包,无监督的离散化可以用infotheo包,有监督的离散化可以用discretization包来处理复杂的离散化操作。

目录

  1. 数据离散化的需求
  2. 无监督的数据离散化
  3. 有监督的数据离散化

1. 数据离散化的需求

数据离散化,教课书上面的定义是,把无限空间中有限的个体映射到有限的空间中。数据离散化操作,大多是针对连续型变量进行的,处理之后数据从连续型变量转变为离散型变量。

离散化处理的特点:

  • 快速迭代,离散化的过程,离散特征增加和减少都很容易,易于扩展。
  • 高效计算,稀疏向量内积乘法运算速度快,计算结果方便存储。
  • 适应性,有些分类模型适合用离散化数据做为数据指标,比如决策树虽然支持输入连续型变量,但决策树模型本身会先将连续型变量转化为离散型变量,做逻辑回归时,也通常会先把连续型特征变量离散化,变成0或1。
  • 鲁棒性,数据离散化会让特征变化不是特别明显,特别是对于降低异常数据干扰。
  • 稳定性,牺牲了数据的微小变化,但保证了模型的稳定性。
  • 过拟合,减少的数据信息,降低了过拟合的风险。
  • 业务性,业务上已经定义出了离散值的含义,需要做离散化处理。

离散化处理通常对连续型变量进行处理,但同时也可以对离散型变量再进行离散化,适用于数据聚合或重新分组。

2. 无监督的数据离散化

我们可以无监督的数据离散化方法,进行数据的分箱操作,包括等宽分箱法、等频分箱法、通过kmeans分箱法等。

  • 等宽分箱法,将观察点均匀划分成n等份,每份的间距相等。
  • 等频分箱法,将观察点均匀分成n等份,每份的观察点数相同。
  • kmeans分箱法,先给定中心数,将观察点利用欧式距离计算与中心点的距离进行归类,再重新计算中心点,直到中心点不再发生变化,以归类的结果做为分箱的结果。

接下来,我们用R语言来实现上面提到的几种分箱的操作。

2.1 准备数据

本文所使用的系统环境

  • Win10 64bit
  • R: 3.2.3 x86_64-w64-mingw32/x64 b4bit

首先,生成一组数据,这组织通过3个正态分布的随机数进行叠加,主要用于体现分箱的不同特征。


# 设置随机数种子,生成符合正态分布N(0,1)的数据1000个点
> set.seed(0)
> a1<-rnorm(1000)

> set.seed(1)
> a2<-rnorm(300,0,0.2)

> set.seed(2)
> a3<-rnorm(300,3,0.5)

# 按顺序,合并3种数据
> aa<-c(a1,a2,a3)

# 查看数据
> head(aa,30)
 [1]  1.26295 -0.32623  1.32980  1.27243  0.41464 -1.53995 -0.92857 -0.29472 -0.00577  2.40465  0.76359
[12] -0.79901 -1.14766 -0.28946 -0.29922 -0.41151  0.25222 -0.89192  0.43568 -1.23754 -0.22427  0.37740
[23]  0.13334  0.80419 -0.05711  0.50361  1.08577 -0.69095 -1.28460  0.04673

画出数据的散点图,x轴为索引序号,y轴为值。


> plot(aa)

查看数据的分布形状,x轴为值,y轴为这个值出现的次数。


# 画出数据的直方图
> hist(aa,200)

通过散点图和直方图,就可以很明显的看到数据的特征,数据有倾斜的,值在0和3附近是比较多的。

2.2 infotheo包使用

在R语言中,我们可以使用infotheo包,来做等宽分箱和等频分箱。

项目主页:https://cran.r-project.org/web/packages/infotheo/

安装infotheo包,是非常简单的,只需要一条命令。


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

进行等宽分箱,将观察点均匀划分成n等份,每份的间距相等。


> d1<-discretize(aa,"equalwidth",5)
> table(d1)
d1
  1   2   3   4   5 
 40 467 703 217 173 

# 可视化分箱
> plot(aa,col=d1$X)

分为5组,第一组最少,第三组最多,每组值的区间是一样的。

进行等频率分箱,将观察点均匀分成n等份,每份的观察点数相同。


> d2<-discretize(aa,"equalfreq",5)
> table(d2)
d2
  1   2   3   4   5 
320 320 320 320 320 

# 可视化分箱
> plot(aa,col=d2$X)

分为5组,每组的数量都相同。

kmeans分箱法,先给定中心数为5。


> d3<-kmeans(aa,5)
> table(d3$cluster)
  1   2   3   4   5 
121 258 303 267 651 

# 5个中心
> d3$centers
     [,1]
1 -1.6196
2  1.2096
3  3.0436
4 -0.7228
5  0.0736

# 可视化分箱
> plot(aa,col=d3$cluster)

5组聚类的结果。

每一种方法,对于分箱的离散化的结果都是不同的,使用哪一种方法,最好都有业务上的解释,不能随便乱用。

3. 有监督的数据离散化方法

有监督的离散化的方法,可以用R语言中discretization包来操作。discretization包,是一个用来做有监督离散化的工具集,主要用于卡方分箱算法,它提供了几种常用的离散化工具函数,可以按照自上而下或自下而上,实施离散化算法。

项目主页:https://cran.r-project.org/web/packages/discretization/

安装discretization包是非常简单的,只需要一条命令。


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

discretization提供了几个主要的离散化的工具函数:

  • chiM,ChiM算法进行离散化
  • chi2, Chi2算法进行离散化
  • mdlp,最小描述长度原理(MDLP)进行离散化
  • modChi2,改进的Chi2方法离散数值属性
  • disc.Topdown,自上而下的离散化
  • extendChi2,扩展Chi2算法离散数值属性

3.1 卡方分箱

卡方分箱是依赖于卡方检验的分箱方法,在统计指标上选择卡方统计量(chi-Square)进行判别。卡方分箱的基本思想是判断相邻的两个区间是否有分布差异,如果两个相邻的区间具有非常类似的分布,则这两个区间可以合并;否则,它们应当保持分开。基于卡方统计量的结果进行自下而上的合并,直到满足分箱的限制条件为止。

卡方分箱的实现步骤:

1. 预先设定一个卡方的阈值。以这个阈值为标准,我们对数据进行卡方检验,通过显著性水平和自由度,计算出数据的卡方值与这个阈值进行比较。

  • 显著性水平,当置信度90%时显著性水平为10%,ChiMerge算法推荐使用置信度为0.90、0.95、0.99。
  • 自由度,比分类数量小1。例如:有3类,自由度为2。

类别和属性独立时,有90%的可能性,计算得到的卡方值会小于4.6。大于阈值4.6的卡方值就说明属性和类不是相互独立的,不能合并。如果阈值选的大,区间合并就会进行很多次,离散后的区间数量少、区间大。

2. 初始化:根据要离散化的数据对实例进行排序,每个实例属于一个区间
3. 合并区间:
1) 计算每一个对相邻区间的卡方值
2) 将卡方值最小的一对区间合并

卡方统计量衡量了区间内样本的频数分布与整体样本的频数分布的差异性,在做分箱处理时可以使用两种限制条件:

  • 分箱个数:限制最终的分箱个数结果,每次将样本中具有最小卡方值的区间与相邻的最小卡方区间进行合并,直到分箱个数达到限制条件为止。
  • 卡方阈值:根据自由度和显著性水平得到对应的卡方阈值,如果分箱的各区间最小卡方值小于卡方阈值,则继续合并,直到最小卡方值超过设定阈值为止。

4.评估指标

分为箱之后需要评估,常用的评估手段是计算出WOE和IV值。对于WOE和IV值的含义,参考文章:https://blog.csdn.net/kevin7658/article/details/50780391

3.2 数据准备

有监督的离散化算法,只少需要数据为2列,一列用于离散化的向量数据,一列是分类的评价标准数据,所以我们需要重新选择一个数据集,具有分类标准的数据。

我们使用一个系统自带的数据集iris,这个着名的Fisher的鸢尾花数据。iris包含150个记录,和5个变量的数据框,名为Sepal.Length(萼片的长度),Sepal.Width(萼片的宽度),Petal.Length(花瓣的长度),Petal.Width(花瓣的长度)和Species(种类,setosa,versicolor和virginica)。

查看数据集


> data(iris)
> head(iris,10)
   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
7           4.6         3.4          1.4         0.3  setosa
8           5.0         3.4          1.5         0.2  setosa
9           4.4         2.9          1.4         0.2  setosa
10          4.9         3.1          1.5         0.1  setosa

接下来,我们分别用不同的算法,把这个数据集进行离散化处理。

3.3 chiM算法进行离散化

ChiM()函数,使用ChiMerge算法基于卡方检验进行自下而上的合并。通过卡方检验判断相邻阈值的相对类频率,是否有明显不同,或者它们是否足够相似,从而合并为一个区间。

chiM(data,alpha)函数解读。

  • 第一个参数data,是输入数据集,要求最后一列是分类属性。
  • 第二个参数alpha,表示显著性水平。
  • 自由度,通过数据计算获得是2,一共3个分类减去1。

下面使用chiM()进行计算。


# 离散化计算
> chi1<-chiM(iris,alpha=0.05)

# 查看前10条结果
> head(chi1$Disc.data,10)
   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1             1           3            1           1  setosa
2             1           2            1           1  setosa
3             1           2            1           1  setosa
4             1           2            1           1  setosa
5             1           3            1           1  setosa
6             1           3            1           1  setosa
7             1           3            1           1  setosa
8             1           3            1           1  setosa
9             1           1            1           1  setosa
10            1           2            1           1  setosa

chi2算法的结果解读,从输出结果可以看出,Sepal.Length Sepal.Width Petal.Length Petal.Width,这几个变量被自动地离散化处理了。

查看这4个列的离散化情况。


> apply(chi1$Disc.data,2,table)
$Sepal.Length
 1  2  3  4 
52 21 65 12 

$Sepal.Width
 1  2  3 
57 56 37 

$Petal.Length
 1  2  3  4 
50 45 21 34 

$Petal.Width
 1  2  3 
50 54 46 

$Species
    setosa versicolor  virginica 
        50         50         50 

再查看每个列的阈值。


# 查看分组阈值
> chi1$cutp
[[1]]
[1] 5.45 5.75 7.05

[[2]]
[1] 2.95 3.35

[[3]]
[1] 2.45 4.75 5.15

[[4]]
[1] 0.80 1.75

以第一组Sepal.Length举例,被分为4类,当源数据Sepal.Length列的值,小于 5.45值为1类,大于5.45同时小于5.75为2类,大于5.75同时小于7.05为3类,大于7.05为4类。

我们把离散化后的数据,与源数据进行合并进行观察。第一列是离散化后的数据,第二列是源数据。


> chi1df<-cbind(chi1$Disc.data$Sepal.Length,iris$Sepal.Length)
> head(chi1df,20)
      [,1] [,2]
 [1,]    1  5.1
 [2,]    1  4.9
 [3,]    1  4.7
 [4,]    1  4.6
 [5,]    1  5.0
 [6,]    1  5.4
 [7,]    1  4.6
 [8,]    1  5.0
 [9,]    1  4.4
[10,]    1  4.9
[11,]    1  5.4
[12,]    1  4.8
[13,]    1  4.8
[14,]    1  4.3
[15,]    3  5.8
[16,]    2  5.7
[17,]    1  5.4
[18,]    1  5.1
[19,]    2  5.7
[20,]    1  5.1

这样就完成了卡方分箱的操作,接下来的其他几个函数的使用,与chiM()的算法类似,就不再过多讨论了。

3.4 其他算法

chi2()算法,查看分组阈值。


> chi2<-chi2(iris,alp=0.5,del=0.05)
> chi2$cutp
[[1]]
[1] 3.5 4.5 6.5

[[2]]
[1] 3.5 4.5

[[3]]
[1] 1.5 2.5 3.5

[[4]]
[1] 1.5 3.5

modChi2()算法,查看分组阈值。


> chi3<-modChi2(iris,alp=0.5)
> chi3$cutp
[[1]]
[1] 1.5 2.5

[[2]]
[1] 2.5

[[3]]
[1] 1.5 2.5

[[4]]
[1] 1.5 2.5

extendChi2()算法,查看分组阈值。


> chi4<-extendChi2(iris,alp = 0.5)
> chi4$cutp
[[1]]
[1] 1.5 2.5

[[2]]
[1] 1.5 2.5

[[3]]
[1] 1.5 2.5

[[4]]
[1] 1.5 2.5

> m1<-mdlp(iris)
> m1$cutp
[[1]]
[1] 5.55 6.15

[[2]]
[1] 2.95 3.35

[[3]]
[1] 2.45 4.75

[[4]]
[1] 0.80 1.75

disc.Topdown()算法,查看分组阈值。


> d1<-disc.Topdown(iris,method=1)
> d1$cutp
[[1]]
[1] 4.30 5.55 6.25 7.90

[[2]]
[1] 2.00 2.95 3.05 4.40

[[3]]
[1] 1.00 2.45 4.75 6.90

[[4]]
[1] 0.10 0.80 1.75 2.50

最后,分箱需要注意的是,分完箱之后,某些箱区间里可能数据分布比例极不均匀,那么这样子会直接导致后续计算WOE时出现inf无穷大的情况,这是不合理的。这种情况,说明分箱太细,需要进一步缩小分箱的数量。

本文详细地介绍了无监督的离散化方法和有监督的离散化方法,针对不同的场景,我们可以选择不同的方法进行使用。R语言中,包提供了各种离散化的工具函数,使用起来很方便,可以大幅提供数据处理过程的效率。

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

打赏作者