• Archive by category "程序算法"

Blog Archives

用R语言10分钟上手神经网络模型

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

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

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

关于作者:

  • 张丹,分析师/程序员/Quant: R,Java,Nodejs
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

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

前言

“神经网络”,“深度学习”,已经不是纯学术的词了,有很多的算法已经落地为我们生活中的各种产品,比如人脸识别,智能客服,从图片中提词,AI作诗,AI作画,chatGPT,AphlaGo等,各种的AI算法应用出圈,打破很多行业壁垒,甚至机器比人类更有创造力,做的更好。这些AI应用出圈的背后,很大程度上是依赖于神经网络的高速发展。

跟上时代的脚步,我将用R语言详细地介绍神经网络建模和应用的过程,以动手操作为主。本文是神经网络的第一篇文章:用R语言10分钟上手神经网络模型neuralnet。

目录

  1. 神经网络介绍
  2. neuralnet包介绍
  3. 神经元模型:零隐藏层
  4. 单层神经网络(感知机):一个隐藏层
  5. 单层神经网络:多分类
  6. 两层神经网络(多层感知器):多分类
  7. 多层神经网络(深度学习)

1. 神经网络介绍

神经网络,也称为人工神经网络 (ANN) 是机器学习的子集,并且是深度学习算法的核心。其名称和结构是受人类大脑的启发,模仿了生物神经元信号相互传递的方式。

人工神经网络 (ANN) 由节点层组成,包含一个输入层、一个或多个隐藏层和一个输出层。 每个节点也称为一个人工神经元,它们连接到另一个节点,具有相关的权重和阈值。 如果任何单个节点的输出高于指定的阈值,那么该节点将被激活,并将数据发送到网络的下一层。 否则,不会将数据传递到网络的下一层。

从单层神经网络(感知器)开始,到包含一个隐藏层的两层神经网络,再到多层的深度神经网络,一共有三次兴起过程。

在神经网络中,处理单元通常按层次分布于神经网络的输入层、隐层和输出层中,因此分别称之为输入节点、隐节点和输出节点,各自的功能如下所示:

  • 输入节点:接受与处理训练数据集中的各输入变量值
  • 隐藏点:实现非线性数据的线性变换
  • 输出节点:给出输出变量的分类或预测结果

2. neuralnet包介绍

在R语言里面,有好几个包都支持神经网络模型,neuralnet包就一个常用包。

neuralnet包,安装过程很简单。


# 安装
> install.packages("neuralnet")

# 加载
> library(neuralnet)

# 设置工作路径
> setwd("C:/work/R/neural/neural_networks")

neuralnet包中的,neuralnet()函数可用于神经网络的模型训练,是这个包中最核心的函数,neuralnet()函数参数很多,是需要我们仔细了解的。

neuralnet()函数,在训练神经网络模型时,支持使用反向传播(BP)、弹性反向传播(RPROP)、无权重回溯、修正的全球收敛版本(GRPROP)训练神经网络,函数允许通过自定义选择误差和激活函数进行灵活设置。此外,还实现了广义权重的计算。

查看neuralnet()函数


> neuralnet
function (formula, data, hidden = 1, threshold = 0.01, stepmax = 1e+05, 
    rep = 1, startweights = NULL, learningrate.limit = NULL, 
    learningrate.factor = list(minus = 0.5, plus = 1.2), learningrate = NULL, 
    lifesign = "none", lifesign.step = 1000, algorithm = "rprop+", 
    err.fct = "sse", act.fct = "logistic", linear.output = TRUE, 
    exclude = NULL, constant.weights = NULL, likelihood = FALSE) 

...

参数列表:

  • formula: 定义模型的公式。
  • data: data.frame格式,训练集原始数据。
  • hidden: 整数向量,指定每层中隐藏神经元的数量,如c(3,2)为二个隐藏层,第一层3个节点,第二层2个节点。
  • threshold: 数值型,停止条件阈值,指定误差函数的偏导数作为停止标准的阈值。
  • stepmax: 数值型,停止条件最大迭代次数,达到这个最大值就会停止神经网络的训练过程。
  • rep: 神经网络训练的重复次数。
  • startweights: 模型起始值的权重向量,设置为NULL用于随机初始化。
  • learningrate.limit: 一个向量或一个列表,包含学习率的最低和最高限制,仅用于RPROP和GRPROP。
  • learningrate.factor: 一个向量或一个列表,包含学习率的上限和下限的乘法系数,仅用于RPROP和GRPROP。
  • learningrate: 数值型,指定传统反向传播使用的学习率,只用于传统的反向传播。
  • lifesign: 字符串,指定函数在计算神经网络时的打印量,’none’, ‘minimal’ or ‘full’。
  • lifesign.step: 整数型,指定在完全lifesign模式下打印最小阈值的步长。
  • algorithm: 字符串,包含用于计算神经网络的算法类型。以下类型是可能的:”backprop”,”rprop+”,”rprop-“,”sag”,或 “slr”。backprop’指的是反向传播,’rprop+’和’rprop-‘指的是带有和不带有权重回溯的弹性反向传播,而’sag’和’slr’则是诱导使用修改后的全局收敛算法(GRPROP)。
  • err.fct: 损失函数,用于计算误差的函数,可以使用字符串’sse’和’ce’,它们代表平方误差之和和交叉熵。
  • act.fct: 激活函数,用于平滑协变量或神经元与权重的交叉积的结果。字符串 “logistic”和 “tanh” 可以用于logistic函数和切线双曲。
  • linear.output:逻辑值,是否线性输出,即是回归还是分类。TRUE表示输出节点的激活函数为线性函数,FALSE表示非线性函数,在B-P中为FALSE
  • exclude: 一个向量或一个矩阵,指定从计算中排除的权重。如果给定的是一个向量,必须知道权重的确切位置。一个有n行3列的矩阵将排除n个权重,其中第一列代表层,第二列代表输入神经元,第三列代表权重的输出神经元。
  • constant.weights: 一个向量,指定训练过程中不需要训练的权重,固定值。
  • likelihood: 逻辑值,如果误差函数等于负对数似然函数,将计算信息准则AIC和BIC。

在实际建模过程中,我们主要使用formula,data,hidden,threshold,algorithm,err.fct,act.fct,linear.output这个参数。

3. 神经元模型:零隐藏层

我们先建立一个最简单的网络模型,只有输入节点和输出节点,不包括隐藏层。这种结构都不算是神经网络,是神经网络的前身神经元模型。神经元模型是一个包含输入,输出与计算功能的模型。输入可以类比为神经元的树突,而输出可以类比为神经元的轴突,计算则可以类比为细胞核。

查看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

以iris数据集的全量数据做为输入,对单一目标Species==”setosa”,做二分类的模型,判断是否属于是setosa的种类。


# 没有隐藏层
> nn <- neuralnet(Species=="setosa"~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width, 
+                 data = iris,hidden = 0)

# 查看网络结构 
> plot(nn)

从输出的神经网络结构图中,我们可以看到,四个输入节点为iris数据集的4个参数(Sepal.Length Sepal.Width Petal.Length Petal.Width Species),输出节点(Species)为模型训练目标,蓝色节点1为这一层的偏差权重(截距)值,从输入节点到输出节点的连线上,有权重值。我们训练神经网络模型,就是为了算出这些权重值,反复迭代。

接下来,我们查看模型的各种输出结果,都在nn这个对象中。

  • call:模型函数
  • response:因变量数据
  • covariate:自变量数据
  • model.list:模型函数
  • err.fct:损失函数
  • act.fct:激活函数
  • linear.output:是否线性输出
  • data:原始数据集
  • exclude:指定从计算中排除的权重
  • net.result:预测值结果
  • weights:各个节点的权重值列表
  • generalized.weights:广义权重
  • startweights:各个节点的初始权重(初始权值为在的正态分布随机数)
  • result.matrix:结果矩阵,终止迭代时各个节点的权重,迭代次数,损失函数值和权重的最大调整量

首先,来查看结果矩阵result.matrix。

> nn$result.matrix
                                             [,1]
error                                1.533438e+00
reached.threshold                    9.212915e-03
steps                                6.297000e+03
Intercept.to.Species == "setosa"     1.202943e-01
Sepal.Length.to.Species == "setosa"  6.432600e-02
Sepal.Width.to.Species == "setosa"   2.441786e-01
Petal.Length.to.Species == "setosa" -2.230368e-01
Petal.Width.to.Species == "setosa"  -5.935682e-02
  • error,损失函数值,1.533438,值越小
  • reached.threshold,终止条件,权重的最大调整量9.212915e-03
  • steps,整个训练执行了6297步
  • Intercept.to.Species,偏差权重(截距),值为1.202943e-01
  • Sepal.Length.to.Species,Sepal.Length变量权重
  • Sepal.Width.to.Species,Sepal.Width变量权重
  • Petal.Length.to.Species,Petal.Length变量权重
  • Petal.Width.to.Species,Petal.Width变量权重

查看连接的初始的权重,随机分配的,分别对应5个不同输入节点到输出节点的权重。


> nn$startweights 
[[1]]
[[1]][[1]]
           [,1]
[1,]  0.5797599
[2,] -0.3526484
[3,]  1.7074446
[4,]  1.3395568
[5,]  0.7781903

查看模型最终的连接权重列表,从startweights初始到weights最终结果。


> nn$weights
[[1]]
[[1]][[1]]
            [,1]
[1,]  0.12029435
[2,]  0.06432600
[3,]  0.24417858
[4,] -0.22303683
[5,] -0.05935682

查看模型中的广义权重,分别对应每个样本数据的权重值


# 取前6条
> head(nn$generalized.weights[[1]])
           [,1]       [,2]       [,3]        [,4]
[1,]  3.1084343  11.799476 -10.777840  -2.8683080
[2,]  0.4883183   1.853634  -1.693141  -0.4505957
[3,]  0.7288879   2.766825  -2.527265  -0.6725814
[4,]  0.4490363   1.704522  -1.556939  -0.4143483
[5,] 20.4485892  77.621919 -70.901170 -18.8689374
[6,] -3.7104548 -14.084718  12.865219   3.4238224

查看模型的计算函数


# 模型公式
> nn$call
neuralnet(formula = Species == "setosa" ~ Sepal.Length + Sepal.Width + 
    Petal.Length + Petal.Width, data = iris, hidden = 0)

# 损失函数
> nn$err.fct
function (x, y) 
{
    1/2 * (y - x)^2
}


attr(,"type")
[1] "sse"

# 激活函数
> nn$act.fct
function (x) 
{
    1/(1 + exp(-x))
}


attr(,"type")
[1] "logistic"

查看模型的变量


# 模型的变量定义
> nn$model.list       
$response
[1] "Species == \"setosa\""

$variables
[1] "Sepal.Length" "Sepal.Width"  "Petal.Length" "Petal.Width" 

# 查看自变量,取前6条
> head(nn$covariate)
     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

查看因变量的实际值,对模型预测输出变量的预测概率值。


# 查看因变量,实际值,取前6条
> head(nn$response)
  Species == "setosa"
1                TRUE
2                TRUE
3                TRUE
4                TRUE
5                TRUE
6                TRUE

# 模型预测结果,取前6条
> head(nn$net.result[[1]])
          [,1]
[1,] 0.9788590
[2,] 0.8439046
[3,] 0.9021788
[4,] 0.8267209
[5,] 0.9968443
[6,] 1.0170459

通过建立一个最简单的网络模型,我们就能了解 neuralnet() 函数的是怎么使用的了,涉及到的内容还是很多的。对于零隐藏层的网络而言,就可以理解为一个线性模型。

4. 单层神经网络(感知机):一个隐藏层

接下来,我们建立一个真正的神经网络模型,包括输入节点、输出节点、1个隐藏层节,还是以iris数据集的全量数据做为输入,对单一目标Species==”setosa”,做二分类的模型,判断是否属于是setosa的种类。

代码只是改动一处,让hidden=1。


# 训练
> n1 <- neuralnet(Species=="setosa"~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width, 
+ data = iris,hidden=1) 

# 画出网络结构 
> plot(n1)

从图的输出,我们看到网络结构中,多了一个隐藏层。

查看结果矩阵


> n1$result.matrix
                                          [,1]
error                             1.899890e-04
reached.threshold                 9.246323e-03
steps                             2.260000e+02
Intercept.to.1layhid1             1.384797e-01
Sepal.Length.to.1layhid1          5.431059e-01
Sepal.Width.to.1layhid1           3.122608e+00
Petal.Length.to.1layhid1         -2.053056e+00
Petal.Width.to.1layhid1          -8.719302e+00
Intercept.to.Species == "setosa" -6.759521e-04
1layhid1.to.Species == "setosa"   1.001522e+00

error的损失值为0.00019,整个训练执行了226步,比神经元模型的零隐藏层,有了大幅的提升。一下子就把分类问题解决了,太神奇了吧。

5. 单层神经网络:多分类

我们增加点难度,用单层神经网路,做个多分类的。


# 单层网络
> n2a <- neuralnet(Species~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width, 
+ data = iris,hidden=1) 

# 画图
> plot(n2a)

从网络结构看,输出为4个节点,隐藏层为一个节点,输出3个节点,为3类。输出结果惨不忍睹,error竟然达到了25,基本都是错了。

接下来,我们试试增加隐藏层节点数,从1个变成2个,看看做个多分类的效果。


# 单层网络,2个节点
> n2b <- neuralnet(Species~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width, 
+ data = iris,hidden=c(2)) 

> plot(n2b)

从网络结构看,输出为4个节点,隐藏层为两个节点,输出3个节点,为3类。输出结果一下就不错了,error为1.9,步骤是15981,这个结果是能用的。

由于一层神经网络的感知器也只能做线性分类任务,在10年后,才对于两层神经网络的研究才带来神经网络的复苏。

6. 两层神经网络(多层感知器):多分类

下面我们要建立一个两层神经网络模型,包括输入层、输出层、2个隐藏层,还是以iris数据集的全量数据做为输入,对单一目标Species做三分类的模型。

首先,先尝2个隐藏层,每层一个节点,试试模型效果。


# 2个隐藏层,每层1个节点
> n3a <- neuralnet(Species~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width, + data = iris,hidden=c(1,1)) > n3a$result.matrix
                                  [,1]
error                     25.000112366
reached.threshold          0.007309118
steps                    358.000000000
Intercept.to.1layhid1      2.677698687
Sepal.Length.to.1layhid1  -0.737245035
Sepal.Width.to.1layhid1    3.844506938
Petal.Length.to.1layhid1  -1.018677154
Petal.Width.to.1layhid1  -10.512582993
Intercept.to.2layhid1     -0.771711595
1layhid1.to.2layhid1       3.469423693
Intercept.to.setosa       -0.510039484
2layhid1.to.setosa         1.612291660
Intercept.to.versicolor    0.754101331
2layhid1.to.versicolor    -0.804506879
Intercept.to.virginica     0.755525924
2layhid1.to.virginica     -0.807147653

error达到了25,又是一个不可接受的结果。

在层数不变的情况下,继续增加每层的节点数,调整为每个隐藏层2个节点。从结果矩阵看,多个很多参数的输出,这些节点都是有神经网络自动生成了,也就是不可解释的因素。


# 2个隐藏层,每层2个节点
> n3b <- neuralnet(Species~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width, + data = iris,hidden=c(2,2)) # 结果矩阵 > n3b$result.matrix
                                 [,1]
error                    24.999473404
reached.threshold         0.008004802
steps                    83.000000000
Intercept.to.1layhid1    -0.493094700
Sepal.Length.to.1layhid1 -1.519304930
Sepal.Width.to.1layhid1  -1.217278013
Petal.Length.to.1layhid1  3.491582287
Petal.Width.to.1layhid1   3.469063900
Intercept.to.1layhid2    -0.921659140
Sepal.Length.to.1layhid2 -2.263419044
Sepal.Width.to.1layhid2   0.410851516
Petal.Length.to.1layhid2  1.581399076
Petal.Width.to.1layhid2  -1.453986728
Intercept.to.2layhid1     0.898644804
1layhid1.to.2layhid1     -2.635824185
1layhid2.to.2layhid1     -2.577338580
Intercept.to.2layhid2    -0.829242612
1layhid1.to.2layhid2     -0.706157301
1layhid2.to.2layhid2      1.549670915
Intercept.to.setosa      -0.436567651
2layhid1.to.setosa        1.527644298
2layhid2.to.setosa        1.165157296
Intercept.to.versicolor   0.374083691
2layhid1.to.versicolor   -1.298521851
2layhid2.to.versicolor    1.805614407
Intercept.to.virginica    0.455763243
2layhid1.to.virginica    -1.181002787
2layhid2.to.virginica     1.253298145

# 画图
> plot(n3b)

效果依然不好,error为24.99,完全不能使用。

下来的操作,就是继续加节点。在层数不变的情况下,继续增加每层的节点数,调整为每个隐藏层3个节点。


# 2个隐藏层,每层3个节点
> n3c <- neuralnet(Species~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width, + data = iris,hidden=c(3,3)) > n3c$result.matrix
                                  [,1]
error                     9.861760e-01
reached.threshold         9.760993e-03
steps                     1.842800e+04
Intercept.to.1layhid1     7.178501e-01
Sepal.Length.to.1layhid1  2.725448e-01
Sepal.Width.to.1layhid1   1.376482e+00
Petal.Length.to.1layhid1 -1.262652e+00
Petal.Width.to.1layhid1  -1.997635e+00
Intercept.to.1layhid2    -2.753156e+00
Sepal.Length.to.1layhid2 -1.465864e+00
Sepal.Width.to.1layhid2   4.125512e+00
Petal.Length.to.1layhid2 -6.026538e-01
Petal.Width.to.1layhid2   2.743376e+00
Intercept.to.1layhid3     1.379399e+00
Sepal.Length.to.1layhid3 -2.256887e-01
Sepal.Width.to.1layhid3   1.027688e+00
Petal.Length.to.1layhid3  3.221391e-01
Petal.Width.to.1layhid3  -2.847147e+00
Intercept.to.2layhid1    -2.353385e+00
1layhid1.to.2layhid1     -1.315389e+03
1layhid2.to.2layhid1      3.555585e+01
1layhid3.to.2layhid1      6.903238e+01
Intercept.to.2layhid2    -6.702997e-01
1layhid1.to.2layhid2      1.352140e+00
1layhid2.to.2layhid2     -5.184099e-02
1layhid3.to.2layhid2     -2.994567e-01
Intercept.to.2layhid3    -3.106167e+00
1layhid1.to.2layhid3      4.958443e+00
1layhid2.to.2layhid3     -2.185125e-01
1layhid3.to.2layhid3     -1.677417e+00
Intercept.to.setosa       4.026726e-01
2layhid1.to.setosa       -1.575905e-02
2layhid2.to.setosa       -1.461912e+00
2layhid3.to.setosa        2.965282e+00
Intercept.to.versicolor   6.790011e-01
2layhid1.to.versicolor   -9.706799e-01
2layhid2.to.versicolor    1.210438e+00
2layhid3.to.versicolor   -2.833834e+00
Intercept.to.virginica   -2.389590e-01
2layhid1.to.virginica     9.950653e-01
2layhid2.to.virginica     7.442383e-01
2layhid3.to.virginica    -3.957608e-01

# 画图
> plot(n3c)

瞬间,模型效果又达到了惊人的新高度,error为0.986,只能用神奇来形容,我也不知道为什么,算法的“玄学” 由此而来。

7. 多层神经网络(深度学习):多个隐藏层

在上个步骤中,我们是层数不变增加每层节点数,当然我们也可以节点数不变,增加更多的层数来操作。继续增加到3层的,保持节点数每个隐藏层2个节点进行测试。


# 3个隐藏层,每层2个节点
> n3d <- neuralnet(Species~ Sepal.Length + Sepal.Width + Petal.Length + 
                         Petal.Width, + data = iris,hidden=c(2,2,2)) 

# 结果矩阵 
> n3d$result.matrix
                                  [,1]
error                     9.814611e-01
reached.threshold         9.856631e-03
steps                     2.599500e+04
Intercept.to.1layhid1     4.377454e+00
Sepal.Length.to.1layhid1  2.988690e+00
Sepal.Width.to.1layhid1   4.895009e+00
Petal.Length.to.1layhid1  4.739920e+00
Petal.Width.to.1layhid1   5.049575e+00
Intercept.to.1layhid2     1.713367e+00
Sepal.Length.to.1layhid2  4.900288e-01
Sepal.Width.to.1layhid2   6.486628e-01
Petal.Length.to.1layhid2 -1.317743e+00
Petal.Width.to.1layhid2  -1.286672e+00
Intercept.to.2layhid1    -5.678148e+01
1layhid1.to.2layhid1     -5.584444e+01
1layhid2.to.2layhid1      9.250502e+02
Intercept.to.2layhid2     1.417653e+01
1layhid1.to.2layhid2      1.370346e+01
1layhid2.to.2layhid2     -3.770254e+01
Intercept.to.3layhid1    -1.702465e+00
2layhid1.to.3layhid1      1.833016e+00
2layhid2.to.3layhid1      4.370202e+00
Intercept.to.3layhid2    -2.525917e+00
2layhid1.to.3layhid2      8.186075e+00
2layhid2.to.3layhid2      1.849959e+00
Intercept.to.setosa       1.991064e+00
3layhid1.to.setosa       -2.193503e+00
3layhid2.to.setosa        1.780463e-01
Intercept.to.versicolor  -2.462871e+00
3layhid1.to.versicolor    2.184454e+00
3layhid2.to.versicolor    1.303726e+00
Intercept.to.virginica    1.471821e+00
3layhid1.to.virginica     9.088568e-03
3layhid2.to.virginica    -1.481810e+00

# 画图
> plot(n3d)

三个隐藏的模型效果,同样能达到非常好的效果,error为0.981。

通过不停的增加层数的节点数,神经网络模型会越来越好,但是计算的复杂度会越来越高,从结果矩阵的输出,我们就可感觉到,输出的变量一次比一次多,都不知道这些变量是干什么用的,要想去解释这个模型根本就无从下手。

保持理性的思维,虽然神经网络可以无限帮我们提升模型的效果,但由于不可以解释,“玄之又玄”,需要在适合场景来使用,而不是全都交给神经网络。本文为神经网络的入门文章第一篇,针对这个主题,我们也准备多写几篇文章,来用R语言进行详细的描述。

说多少理论都不如上手一试。本文代码已上传到github: https://github.com/bsspirit/neural_networks/blob/main/01_neuralnet.r

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

打赏作者

知识图谱:最短路径算法Bellman-Ford

算法为王系列文章,涵盖了计算机算法,数据挖掘(机器学习)算法,统计算法,金融算法等的多种跨学科算法组合。在大数据时代的背景下,算法已经成为了金字塔顶的明星。一个好的算法可以创造一个伟大帝国,就像Google。

算法为王的时代正式到来….

关于作者:

  • 张丹, 分析师/程序员/Quant: R,Java,Nodejs
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-graph-shortest-path-bellman-ford/

前言

最短路径是图算法中的基本算法,被普遍应用于多种场景中,如企业生产经营的场景中,现金流的生产和交换的过程,会产生正向的现金流比如经营收入所得,也会产生负向的现金流比如借款、支出等。如果我们以企业作为节点,以现金流做为边时,由于出现了负向的现金流的边,这样构成负权重的网络图谱,就需要用Bellman-Ford算法来实现。

常用的最短路径算法包括:Dijkstra算法,A star算法,Bellman-Ford算法,SPFA算法(Bellman-Ford算法的改进版本),Floyd-Warshall算法,Johnson算法以及Bi-direction BFS算法。

目录

  1. Bellman-Ford算法介绍
  2. R语言调用Bellman-Ford算法
  3. 自己手写Bellman-Ford算法
  4. 算一个复杂点的网络关系图

1. Bellman-Ford算法介绍

Bellman Ford算法最初于1955年被Alfonso Shimbel提出,但最终基于1958和1956年发表论文的 Richard Bellman和Lester Ford, Jr.二人命名。Edward F. Moore于1957年也提出了该算法,因此该算法有时也称作是Bellman Ford Moore算法。

Bellman-Ford算法与Dijkstra算法一样,解决的是单源最短路径问题。两者不同之处在于,Dijkstra只适用于无负权边的图,而Bellman-Ford算法无此限制:无论是否有向,无论是否有负权边,只要图中没有负权环,则该算法可以正确地给出起点到其余各点的最短路径,否则报告负权环的存在。

Bellman-Ford 算法采用动态规划进行设计,实现的时间复杂度为 O(V*E),其中 V 为顶点数量,E 为边的数量。每次松弛操作实际上是对相邻节点的访问(相当于广度优先搜索),第n次松弛操作保证了所有深度为n的路径最短。由于图的最短路径最长不会经过超过|V| – 1条边,所以可知贝尔曼-福特算法所得为最短路径,也可只时间复杂度为O(VE)。

Bellman-Ford 算法描述:

  • 创建源顶点 v 到图中所有顶点的距离的集合 distSet,为图中的所有顶点指定一个距离值,初始均为 Infinite,源顶点距离为 0;
  • 计算最短路径,执行 V – 1 次遍历;
  • 对于图中的每条边:如果起点 u 的距离 d 加上边的权值 w 小于终点 v 的距离 d,则更新终点 v 的距离值 d;
  • 检测图中是否有负权边形成了环,遍历图中的所有边,计算 u 至 v 的距离,如果对于 v 存在更小的距离,则说明存在环;

下面,我们来举例说明一下算法,创建一个带负权重的有向图,来描述对最短路径

开始点为A,数组U为开始列,数组S为上层节点列。通过6次遍历,找到从A到G的最短路径。

通过4次遍历,其实就能获得从A到D的最短路径。

2. R语言调用Bellman-Ford算法

接下来,我们可以继续使用R语言的igraph包来进行实现。新建数据集,如上面例子所示。

# 加载igraph库
> library(igraph)

# 新建图数据
> df<-data.frame(
+   from=c("A","A","B","B","B","D","D","E"),
+   to=c("B","C","C","D","E","B","C","D"),
+   weight=c(-1,4,3,2,2,1,5,-3)
+ )
 
# 用igraph画图 
> a<-graph_from_data_frame(df, directed = TRUE)
> E(a)$arrow.size<-1
> E(a)$arrow.width<-0.2
> E(a)$label<-df$weight

# layout需要设置一下,不然图的展示会很不好看
> plot(a,layout=layout.gem)

layout需要设置一下,不然图的展示会节点叠加在一起,就看不出来图的样子了,我这里选择layout.gem类型。

然后,使用shortest.paths()函数,进行bellman-ford算法计算最短路径的时候,就报错了。错误提示,图形成负权重的环路。

# 使用bellman-ford算法
> shortest.paths(a,"A", algorithm = "bellman-ford")
Error in shortest.paths(a, "A", algorithm = "bellman-ford") : 
  At core/paths/bellman_ford.c:157 : cannot run Bellman-Ford algorithm, Negative loop detected while calculating shortest paths

尝试了多种数据集,只要是有负值shortest.paths()函数都会报错,无论是否有负权重的环路,我也尝试找了几个的其他的R包,竟然没有发现Bellman-Ford算法的实现,因此只好自己手动写一个了。

3. 自已手写Bellman-Ford算法

于是,我就自己实现bellman-ford算法函数。写了2个小时,好像不练习了,可能还是会有bug吧。

在实现过程中,我们没有做“环”的判断,可能会导致计算逻辑出现问题。我假设了节点的寻路的过程是排序好的,不会出现对于一节点的反复寻路,当然这个假设条件,在某些复杂的场景中是不靠谱的,所以大家在使用的我的代码的时候,当一个基本版参考和学习用,不要直接用于自己的生产环境。

# bellman-ford 最短距离算法
# @param df数据集
# @param start开始节点
bellman_ford<-function(df,start="A"){
  
  # 初始化结果数据集
  rs<-data.frame(from=start,to=start,weight=0)
  
  # 初始化结束节点
  finishnode<-c()
  
  # 初始化开始节点
  nodes<-rs
  
  # 当开始节点数量<结束节点数量时,继续循环计算
  while(length(finishnode)<length(rs$to)){
    
    # 过滤已计算过的结束节点
    if(length(finishnode)>0){
      nodes<-rs[-which(rs$to %in% finishnode),]  
    }
    
    # 对原始数据进行循环,遍历每一行
    for(i in 1:nrow(df)){
      row<-df[i,];row
      
      # 对结束节点进行循环,遍历每一个结束节点
      for(j in nodes$to){
        if(row$from==j){
          # 把之前结束节点的权重叠加本次计算中
          row$weight<-row$weight+rs$weight[which(rs$to==j)]
          # 把结束节点加入到结果数据集
          rs<-rbind(rs,row)
        } 
      }
    }
    
    # 保留最短的唯一结束节点
    rs2<-rs[order(rs$weight),]
    n<-which(duplicated(rs2$to))
    if(length(n)<0){
      rs2<-rs2[-n,]
      rs<-rs2[order(as.numeric(row.names(rs2))),]
    }
    
    # 更新已经完成的结束节点
    finishnode<-c(finishnode,nodes$to)
  }
  return(rs)
}

执行计算,计算A点其他所有节点的距离。

> rs<-bellman_ford(df,"A")
> rs
  from to weight
1    A  A      0
2    A  B     -1
3    B  C      2
5    B  E      1
8    E  D     -2

结果数据集中,to 表示从开始点A到节点的距离,每一行from 到 to 表示,最短路径的顺序。

# 计算2点间的最短路径
# @param rs 结果数据集
# @param start 开始节点
# @param end 结束节点
get_shortest_path<-function(rs,start,end){
  
  # 初始化节点
  q<-c(start)
  idx<-which(rs$from==start)
  
  # 判断有几个分支
  if(length(idx)>0){
    for(i in idx){
      # 递归计算
      q2<-get_shortest_path(rs,rs$to[i],end)
      if(tail(q2,1)==end){
        q<-c(q,q2)
      }
    }
  }
  
  # 查看队列状态
  # print(q)
  
  # 当找到结束节点,程序停止
  if(tail(q,1)==end){
    return(q)
  }
  return(q)
}

# 计算从开始节点到结束节点的路径
> path<-get_shortest_path(rs[-1,],"A","D")
> path
[1] "A" "B" "E" "D"

最后,我们把从A点到D的最短路径涂上颜色,来出来。

# 把路径合并成data.frame
> path2<-data.frame(from=path[-length(path)],to=path[-1])
# 找到最短路径的索引
> idx<-which(paste0(df$from,"-",df$to) %in% paste0(path2$from,"-",path2$to))

# 涂色
> df$color<-"gray"
> df$color[idx]<-"red"
> E(a)$color<-df$color
> plot(a,layout=layout.gem)

红色路径就表示从A到D点的最短路径,看上去结果也是对的。在实现过程中,我们没有做“环”的判断,可能会导致计算逻辑出现问题。自己用写了一段算法,用R语言写算法,感觉比起C或者JAVA还要是简单很多的,因为大量基础函数提供了很多的功能和数据结构支持,不用像底层语言什么结构都需要自己从头写起。

4. 算一个复杂点的网路关系图

生成新数据集,共11个节点。

> df2<-data.frame(
+   from=c("A","B","C","D","C","F","B","B","F","G","H","F","E"),
+   to=  c("B","C","D","C","E","E","F","G","G","H","I","I","J"),
+   weight=c(5,20,10,-5,75,25,30,60,5,-50,-10,50,100)
+ )

用igraph进行画图展示。

> a2<-graph_from_data_frame(df2, directed = TRUE)
> E(a2)$arrow.size<-1
> E(a2)$arrow.width<-0.2
> E(a2)$label<-df2$weight
> plot(a2,layout=layout.gem)

计算A开始到所有节点的最短路径。

> rs<-bellman_ford(df2,"A")
> rs
   from to weight
1     A  A      0
2     A  B      5
3     C  D     35
6     F  E     60
7     B  F     35
9     F  G     40
10    G  H    -10
11    H  I    -20
13    E  J    160
21    B  C     25

计算A到D的最短路径路线,并画科展示。

# 计算 A 到 D的最短路径
> path<-get_shortest_path(rs[-1,],"A","D")
> path
[1] "A" "B" "C" "D"

# 画图
> path2<-data.frame(from=path[-length(path)],to=path[-1])
> idx<-which(paste0(df2$from,"-",df2$to) %in% paste0(path2$from,"-",path2$to))
> df2$color<-"gray"
> df2$color[idx]<-"red"
> E(a2)$arrow.size<-1
> E(a2)$arrow.width<-0.2
> E(a2)$color<-df2$color
> plot(a2,layout=layout.gem)

计算A到I的最短路径路线,并画科展示。

# 计算 A 到 I的最短路径
> path<-get_shortest_path(rs[-1,],"A","I")
> path
[1] "A" "B" "F" "G" "H" "I"

# 画图
> path2<-data.frame(from=path[-length(path)],to=path[-1])
> idx<-which(paste0(df2$from,"-",df2$to) %in% paste0(path2$from,"-",path2$to))
> df2$color<-"gray"
> df2$color[idx]<-"red"
> E(a2)$arrow.size<-1
> E(a2)$arrow.width<-0.2
> E(a2)$color<-df2$color
> plot(a2,layout=layout.gem)

算法和数据结构还是很重要的,图论算法原来仅在计算领域有应用,随着知识图谱的应用普及,相信未来会有更广阔的应用前景,把原来藏在背后的知识,扩展到前台形成新的模型思路。

努力吧,曾经的少年。

本文算法代码已上传到github: https://github.com/bsspirit/graph/blob/main/Bellman-Ford.r

转载请注明出处:
http://blog.fens.me/r-graph-shortest-path-bellman-ford/
打赏作者

知识图谱:最短路径算法Dijkstra

算法为王系列文章,涵盖了计算机算法,数据挖掘(机器学习)算法,统计算法,金融算法等的多种跨学科算法组合。在大数据时代的背景下,算法已经成为了金字塔顶的明星。一个好的算法可以创造一个伟大帝国,就像Google。

算法为王的时代正式到来….

关于作者:

  • 张丹, 分析师/程序员/Quant: R,Java,Nodejs
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-graph-shortest-path-dijkstra/

前言

最短路径是图算法中的基本算法,被普遍应用于多种场景中,比如在实际生活中的路径规划、地图导航,建立道路交通网、物流运输网络、计算机网络等,这时就可以使用最短路径算法。

常用的最短路径算法包括:Dijkstra算法,A star算法,Bellman-Ford算法,SPFA算法(Bellman-Ford算法的改进版本),Floyd-Warshall算法,Johnson算法以及Bi-direction BFS算法。

目录

  1. 最短路径算法介绍
  2. Dijkstra 算法原理
  3. 用R语言调用 Dijkstra 算法

1. 最短路径算法介绍

最短路径问题是图论研究中一个经典算法问题,旨在寻找图中两节点或单个节点到其他节点之间的最短路径。根据问题的不同,算法的具体形式包括:

  • 确定起点的最短路径问题,即给定起始节点,求该节点到其他剩余节点的最短路径,适合使用Dijkstra算法
  • 确定终点的最短路径问题,即给定终点,求其他节点到该终点的最短路径。在无向图中,该问题与确定起点的问题等价;在有向图中,问题等价于把所有路径方向反转的确定起点的问题;
  • 确定起点终点的最短路径问题,即给定起点和终点,求两节点之间的最短路径;
  • 全局最短路径问题,即求图中所有节点之间的最短路径,适合使用Floyd-Warshall算法。

多种最短路径算法

最短路径算法描 述
迪杰斯特拉算法
Dijkstra
寻找某个特定顶点到其它所有顶点的最短路径,该算法要求所有路径的权值为非负数。
贝尔曼福特算法
Bellman-Ford
寻找某个特定顶点到其它所有顶点的最短路径,该算法允许路径的权值为负数。
弗洛伊德算法
(Floyd-Warshall)
寻找各个顶点之间的最短路径,允许非环路的路径权值为负数,该算法不仅适用于稀疏图,在稠密图(路径数量多的图)中寻找最短路径的效率也很高。
约翰逊算法
(Johnson)
寻找各个顶点之间的最短路径,允许非环路的路径权值为负数,该算法更适用于稀疏图(路径数量少的图)。

2. Dijkstra 算法原理

Dijkstra算法,于1956年由荷兰计算机科学家艾兹赫尔.戴克斯特拉提出,用于解决赋权有向图的单源最短路径问题。所谓单源最短路径问题是指确定起点,寻找该节点到图中任意节点的最短路径,算法可用于寻找两个城市中的最短路径或是解决著名的旅行商问题。

基本思想

  1. 通过Dijkstra计算图中的最短路径时,需要指定一个起点A。
  2. 简历两个数组S和U,S用于记录已求出最短路径的顶点(以及相应的最短路径长度),U用于记录还未求出最短路径的顶点(以及该顶点到起点A的距离)。
  3. 初始时,数组S中只有起点A;数组U中是除起点A之外的顶点,并且数组U中记录各顶点到起点A的距离。如果顶点与起点A不相邻,距离为无穷大INF。
  4. 从数组U中找出路径最短的点X,并将其加入到数组S中;同时,从数组U中移除顶点X。接着,更新数组U中的各顶点到起点A的距离。
  5. 重复第4步操作,直到遍历完所有顶点。

下面,我们来举例说明一下算法,创建一个带权重的有向图,来描述对最短路径。

开始点为A,数组U为开始列,数组S为上层节点列。通过6次遍历,找到从A到G的最短路径。

3. 用R语言调用 Dijkstra 算法

我们用R语言中,igraph包已经实现了Dijkstra算法,我们可以直接调用igraph包的shortest.pahs()函数进行最短路径的计算。首先,安装和加载 igraph 程序包

> install.pacakges("igraph")
> library(igraph)

生一个有向有权重的网络图,如上例所示。

> df<-data.frame(
+   from=c("A","A","A","B","D","B","C","D","C","F","E","F"),
+   to=c("B","D","C","C","C","E","E","F","F","E","G","G"),
+   weight=c(4,6,6,1,2,7,6,5,4,1,6,8)
+ )
> df
   from to weight
1     A  B      4
2     A  D      6
3     A  C      6
4     B  C      1
5     D  C      2
6     B  E      7
7     C  E      6
8     D  F      5
9     C  F      4
10    F  E      1
11    E  G      6
12    F  G      8

把数据生成igraph图对象,再以网络图进行展示。

# 生成igraph对象
> a<-graph_from_data_frame(df, directed = TRUE, vertices = NULL)

# 画图
> E(a)$arrow.size<-1
> E(a)$arrow.width<-0.2
> E(a)$label<-df$weight
> plot(a)

上面我们把节点涂色了,下面继续把边也涂色。首先,找到最短路径上的边,把最短路径边涂红色,非最短路径的边涂灰色。

# 通过最短路径的节点,建立边的关系
> s_vert<-data.frame(from=s_node[-length(s_node)],to=s_node[-1])

# 找到网络图中所有的边
> a_vert<-as_data_frame(a,what="edges")

# 找到最短路径的边的索引
> idx<-which(paste0(a_vert$from,"-",a_vert$to) %in% paste0(s_vert$from,"-",s_vert$to))
> idx
[1]  1  4  9 10 11

# 涂色:红色为最短路径的边,灰色为非最短路径的边
> a_vert$color<-"gray"
> a_vert$color[idx]<-"red"

# 画图
> E(a)$color<-a_vert$color
> plot(a)

用Dijkstra算法,计算从A开始,到B,C,D,E,F,G的所有节点的最短距离。

> shortest.paths(a,"A", algorithm = "dijkstra")
  A B D C F  E  G
A 0 4 6 5 9 10 16

输出结果,从A到B最短距离为4,从A到D最短距离为6,从A到C最短距离为5,从A到F最短距离为9,从A到E最短距离为10,从A到G最短距离为16。

当然,我们除了想知道最短距离,还想知道最短距离的路径是什么样的。接下来,我们先计算最短距离路径上的点,并发点进行涂色。

# 计算最短距离经过的点
> s_node<-all_shortest_paths(a,"A","G")$res[[1]]$name
> snode
[1] "A" "B" "C" "F" "E" "G"

# 把节点进行涂色
> cols<-data.frame(label=V(a)$name,color=0)
> cols$color[which(cols$label %in% s_node)]<-1
> cols
  label color
1     A     1
2     B     1
3     D     0
4     C     1
5     F     1
6     E     1
7     G     1

# 画图
> V(a)$color<-cols$color
> plot(a)

其中,黄色点为最短距离穿过的节点,白色节点为不在最短距离路径上的节点。

上面我们把点进行了涂色,下面我们继续边最少路径的边进行涂色。首先,根据最短路径的节点,创建边的关系。然后对边进行涂色,在最短路径的边涂为红色,不在最短路径的边涂为灰色。

# 建立最短路径节点的关系
> s_vert<-data.frame(from=s_node[-length(s_node)],to=s_node[-1])

# 找到网络图的完整路径
> a_vert<-as_data_frame(a,what="edges")

# 找到最短路径的边的索引
> idx<-which(paste0(a_vert$from,"-",a_vert$to) %in% paste0(s_vert$from,"-",s_vert$to))
> idx
[1]  1  4  9 10 11

# 涂色
> a_vert$color<-"gray"
> a_vert$color[idx]<-"red"

# 画图
> E(a)$color<-a_vert$color
> plot(a)

这样就在完整的网络图中,可以从开始节点A,按红色路径走到G,生成了最短路径所过过的节点和边了。

扩展一下,我们用Dijkstra算法,分别计算从不同点开始,到所有不同的点的最短距离,通过矩阵来表达,任意2个节点的最短距离。

> shortest.paths(a, algorithm = "dijkstra")
   A  B  D  C F  E  G
A  0  4  6  5 9 10 16
B  4  0  3  1 5  6 12
D  6  3  0  2 5  6 12
C  5  1  2  0 4  5 11
F  9  5  5  4 0  1  7
E 10  6  6  5 1  0  6
G 16 12 12 11 7  6  0

通过上面的步骤,我们就实现了在R语言中使用Dijkstra算法,计算的最短路径算法。本来还想着,自己再来实现一遍,发现网上有很多的开源的代码,还让自己休息一下吧。下篇文章将介绍,知识图谱:最短路径算法Bellman-Ford算法,应该是需求自己手写了。

本文算法代码已上传到github: https://github.com/bsspirit/graph/blob/main/dijkstra.r

转载请注明出处:
http://blog.fens.me/r-graph-shortest-path-dijkstra/
打赏作者

线性和非线性的最小二乘回归

算法为王系列文章,涵盖了计算机算法,数据挖掘(机器学习)算法,统计算法,金融算法等的多种跨学科算法组合。在大数据时代的背景下,算法已经成为了金字塔顶的明星。一个好的算法可以创造一个伟大帝国,就像Google。

算法为王的时代正式到来….

关于作者:

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

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

前言

很多经典的机器学习模型都是在最小二乘法的基础上发展起来的,理解最小二乘法的原理,对于我们掌握机器学习的高级算法是非常重要的,尤其是回归模型中。最小二乘法是主要是用来做函数拟合,或者求函数极值等问题。

本文以回归模型为例,介绍在最小二乘法在线性和非线性模型中的应用。

目录

  1. 最小二乘法介绍
  2. 最小二乘求解
  3. 线性最小二乘回归
  4. 非线性最小二乘回归

1. 最小二乘法介绍

最小二乘法(Least Squares Method,简记为LSE),源于天文学和测地学上的应用需要,是勒让德( A. M. Legendre)于1805年在其著作《计算慧星轨道的新方法》中提出的,它的主要思想是通过计算理论值与观测值之差的平方和最小,求解未知的参数。

最小二乘法的公式为:

公式解读:

  • H, 目标值,观测值与理论值之差的平方和
  • H’, 最小化的H值
  • y, 观测值
  • yi, 理论值
  • argmin, 最小化函数

我们先了把术语进行统一定义,观测值通常说的是我们的样本数据,理论值就是假设拟合函数产生的预测数据,目标函数是判断一个模型的好不好的标准,进行最小二乘法建立目标函数,让目标函数最小化。

以上图为线性拟合为例,观测值为y1,y2,…,y7,理论值为y1′,y2′,….,y7’,是线上的点并与观测值的x相同,红色线就是我们要做的拟合直线,我们可以画出无数条这样的直线,判断标准是让这条线上的点的理论值与观测值的误差平方和最小,即最小二乘法。

SSE=sum((y1-y1')^2+(y2-y2')^2+(y3-y3')^2+(y4-y4')^2+(y5-y5')^2+(y6-y6')^2+(y7-y7')^2)

利用最小二乘法,可以对数据进行线性和非线性拟合,使误差平方和(SSE)或残差平方和最小。如果观测到的误差近似正态分布时,这种方法是非常有效的。我们在线性回归模型中,要进行残差的正态分布检验,就是基于这个目的的,可以参考文章R语言解读一元线性回归模型

2. 最小二乘求解

在线性方程的求解或是数据曲线拟合中,利用最小二乘法求得的解则被称为最小二乘解。最小二乘法(又称最小平方法)是一种数学优化技术,它通过最小化误差的平方和寻找数据的最佳函数匹配。其他一些优化问题也可通过最小化能量或最大化熵用最小二乘法来表达。

定义数据集x和y,两个变量。

> x<-c(6.19,2.51,7.29,7.01,5.7,2.66,3.98,2.5,9.1,4.2)
> y<-c(5.25,2.83,6.41,6.71,5.1,4.23,5.05,1.98,10.5,6.3)

可视化,画出x,y变量的散点图。

> plot(x,y)

对数据进行最小二乘求解,找到最小二乘解。


> fit<-lsfit(x,y);fit
$coefficients
Intercept         X
0.8310557 0.9004584

$residuals
[1] -1.1548933 -0.2612063 -0.9853975 -0.4332692 -0.8636686  1.0037250  0.6351198 -1.1022017
[9]  1.4747728  1.6870190

$intercept
[1] TRUE

$qr
$qt
[1] -17.1901414   6.2044421  -0.7047339  -0.1530724  -0.5856560   1.2766692   0.9102648
[8]  -0.8295242   1.7584540   1.9625308

$qr
Intercept           X
[1,] -3.1622777 -16.1718880
[2,]  0.3162278   6.8903149
[3,]  0.3162278  -0.2782874
[4,]  0.3162278  -0.2376506
[5,]  0.3162278  -0.0475287
[6,]  0.3162278   0.3936703
[7,]  0.3162278   0.2020970
[8,]  0.3162278   0.4168913
[9,]  0.3162278  -0.5409749
[10,]  0.3162278   0.1701682

$qraux
[1] 1.316228 1.415440

$rank
[1] 2

$pivot
[1] 1 2

$tol
[1] 1e-07

attr(,"class")
[1] "qr"

结果解读:

  • coefficients, 系数的最小二乘估计,包括截距Intercept和自变量X
  • residuals, 残差
  • intercept, 有截距
  • qt, 矩阵QR分解。QR分解法是将矩阵分解成一个正规正交矩阵与上三角形矩阵。关于矩阵的详细操作,请参考文章R语言中的矩阵计算

我们看到lsfit()函数使用矩阵的QR分解的方法,计算得出了最小二乘估计的截距Intercept和自变量系数。

目标:y=a*x+b,用最小二乘法估计出参数a和b。

根据最小二乘法的矩阵计算公式,我直接跳过了推到过程,得出结论。

公式解释:

  • X,为自变量矩阵
  • Y,为因变量矩阵
  • B,为参数,需要求解
  • XT,X的转置矩阵
  • *, 矩阵乘法
  • -1, 计算逆矩阵

手动计算过程:


# 因变量矩阵
> y1<-matrix(y,ncol=1)  

# 自变量矩阵,增加截距列
> x1<-matrix(c(x,rep(1,length(x))),ncol=2)

# 求解最小二乘解
> solve(t(x1)%*%x1) %*% t(x1) %*% y1
          [,1]
[1,] 0.9004584
[2,] 0.8310557

对应到一元线性回归方程中,y = 0.9004584 * x + 0.8310557,与lsfit()的结果一致。

3. 线性最小二乘回归

接下来,我们先进行线性回归分析,以上面的数据集作为样本,只包含x和y。以普通最小二乘法进行线性回归,计算出预测值并和观测值进行比较,计算预测值和观测值的误差平方和,使得误差函数最小。线性回归模型的详细介绍,可以参考文章R语言解读一元线性回归模型, R语言解读多元线性回归模型

建立线性回归模型


# 建立回归模型
> line<-lm(y~x)

# 查看回归模型
> summary(line)
Call:
lm(formula = y ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.1549 -0.9550 -0.3472  0.9116  1.6870 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.8311     0.9440   0.880 0.404338    
x             0.9005     0.1698   5.302 0.000726 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.17 on 8 degrees of freedom
Multiple R-squared:  0.7785,	Adjusted R-squared:  0.7508 
F-statistic: 28.12 on 1 and 8 DF,  p-value: 0.0007263

# 残差平方和
> sum(line$residuals^2)
[1] 10.95334

通过查看模型的结果数据,我们可以发现通过T检验自变量x是非常显著,但截距Intercept不显著。通过F检验判断出整个模型的自变量方差是非常显著。通过R^2的相关系数检验可以判断自变量和因变量是不是高度相关的,残差平方和(residual sum-of-squares)为10.95334,效果不好。

接下来,使用模型进行预测,生成预测值。


# 设置x的取值区间
> new <- data.frame(x = seq(min(x),max(x),len = 100))

# 预测结果
> pred<-predict(line,newdata=new)

进行可视化,把观测值与预测值进行展示,以散点表现观测值,以线表示预测值。


> plot(x,y)
> lines(new$x,pred,col="red",lwd=2)

从图中我们可以很容易看出,拟合的效果并不好,与R^2检验不通过是一致的。那么,如果线性模型表现不好,我们可以试试非线性的模型,还是用最小二乘法进行非线性的回归。

4. 非线性最小二乘回归

由于散点的分布并不是线性的,我们再尝试用非线性的函数,进行最小二乘回归的拟合,包括一元二次函数,一元三次函数,指数函数,倒数函数。在R语言中,我们可以使用nls()函数,进行非线性最小二乘回归。

4.1 一元二次函数非线性拟合
定义公式 y ~ a*x+b*x^2+c ,其中x, x^2是一次函数和二次函数,a,b,c分别是参数,我们需要给定参数的初始值。


# 建立回归模型
> nline1 <- nls(y ~ a*x+b*x^2+c, start = list(a = 1,b = 1,c=2))

# 查看模型结果
> nline1
Nonlinear regression model
  model: y ~ a * x + b * x^2 + c
   data: parent.frame()
       a        b        c 
0.008987 0.082188 2.850389 
 residual sum-of-squares: 9.758

Number of iterations to convergence: 1 
Achieved convergence tolerance: 3.024e-06

# 分析模型
> summary(nline1)

Formula: y ~ a * x + b * x^2 + c

Parameters:
  Estimate Std. Error t value Pr(>|t|)
a 0.008987   0.977890   0.009    0.993
b 0.082188   0.088760   0.926    0.385
c 2.850389   2.379763   1.198    0.270

Residual standard error: 1.181 on 7 degrees of freedom

Number of iterations to convergence: 1 
Achieved convergence tolerance: 3.024e-06

通过查看模型的结果数据,我们可以发现通过T检验相关系数a,b,c都不显著,残差平方和(residual sum-of-squares)为9.758,数值过大,效果不好。

把观测值与预测值进行可视化展示,以散点表现观测值,以线表示预测值。


# 模型预测
> npred1<-predict(nline1,newdata=new)

# 可视化
> plot(x,y)
> lines(new$x,npred1,col="blue",lwd=2)

从图中我们可以看出,拟合的是一条曲线,效果不太好。接下来,我们再试试一元三次函数。

4.2 一元三次函数非线性拟合

定义公式 y ~ a*x+b*x^2+c*x^3+d ,其中x, x^2, x^3是一次函数,二次函数和三次函数,a,b,c,d分别是参数,我们需要给定参数的初始值。

使用一元三次函数,进行非线性拟合。


> nline2 <- nls(y ~ a*x+b*x^2+c*x^3+d, start = list(a = 1,b = 1,c=2,d=1))

# 查看模型结果
> nline2
Nonlinear regression model
  model: y ~ a * x + b * x^2 + c * x^3 + d
   data: parent.frame()
      a       b       c       d 
 10.011  -1.822   0.110 -12.516 
 residual sum-of-squares: 3.453

Number of iterations to convergence: 2 
Achieved convergence tolerance: 3.983e-08

# 分析模型
> summary(nline2)

Formula: y ~ a * x + b * x^2 + c * x^3 + d

Parameters:
   Estimate Std. Error t value Pr(>|t|)  
a  10.01051    3.08630   3.244   0.0176 *
b  -1.82152    0.57797  -3.152   0.0198 *
c   0.10998    0.03323   3.310   0.0162 *
d -12.51586    4.88778  -2.561   0.0429 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7586 on 6 degrees of freedom

Number of iterations to convergence: 2 
Achieved convergence tolerance: 3.983e-08

通过查看模型的结果数据,我们可以发现通过T检验相关系数a,b,c,d在0.01的置信区间显著,残差平方和(residual sum-of-squares)为3.453,数值较小,效果应该算是还可以的,能够表达拟合的特征。

把观测值与预测值进行可视化展示,以散点表现观测值,以线表示预测值。


> npred2<-predict(nline2,newdata=new)
> plot(x,y)
> lines(new$x,npred2,col="green",lwd=2)

从图中我们可以看出,拟合的曲线基本穿过观测值的区域,直观上说明效果还不错,与数据检验给出的结论是一致的。

4.3 指数函数非线性拟合

定义公式 y ~ a*exp(b*x)+c ,其中exp(x)是指数函数,a,b,c分别是参数,我们需要给定参数的初始值。


> nline3 <- nls(y ~ a*exp(b*x)+c, start = list(a = 1,b = 1,c=2))

> nline3
Nonlinear regression model
  model: y ~ a * exp(b * x) + c
   data: parent.frame()
     a      b      c 
0.6720 0.2718 2.2087 
 residual sum-of-squares: 8.966

Number of iterations to convergence: 17 
Achieved convergence tolerance: 7.136e-06

> summary(nline3)

Formula: y ~ a * exp(b * x) + c

Parameters:
  Estimate Std. Error t value Pr(>|t|)
a   0.6720     1.2127   0.554    0.597
b   0.2718     0.1764   1.541    0.167
c   2.2087     2.2447   0.984    0.358

Residual standard error: 1.132 on 7 degrees of freedom

Number of iterations to convergence: 17 
Achieved convergence tolerance: 7.136e-06

通过查看模型的结果数据,我们可以发现通过T检验相关系数a,b,c都不显著,残差平方和(residual sum-of-squares)为8.966,数值过大,效果不好。


> npred3<-predict(nline3,newdata=new)
> plot(x,y)
> lines(new$x,npred3,col="yellow",lwd=2)

从图中我们可以看出,指数函数的拟合曲线,和二次函数表现差不多,效果不太好。

4.4 倒数函数非线性拟合

定义公式 y ~ y ~ b/x+c ,其中1/x是倒数函数,b,c分别是参数,我们需要给定参数的初始值。

使用倒数函数,进行拟合。

> nline4 <- nls(y ~ b/x+c, start = list(b = 1,c=2))

# 查看模型结果
> nline4
Nonlinear regression model
  model: y ~ b/x + c
   data: parent.frame()
    b     c 
-17.0   9.5 
 residual sum-of-squares: 15.74

Number of iterations to convergence: 1 
Achieved convergence tolerance: 1.259e-07

> summary(nline4)

Formula: y ~ b/x + c

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
b  -17.003      4.108  -4.139  0.00326 ** 
c    9.500      1.077   8.817 2.15e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.403 on 8 degrees of freedom

Number of iterations to convergence: 1 
Achieved convergence tolerance: 1.259e-07

通过查看模型的结果数据,我们可以发现通过T检验相关系数b,c非常显著,残差平方和(residual sum-of-squares)为15.74,数值较大,效果不好。

把观测值与预测值进行可视化展示,以散点表现观测值,以线表示预测值。


> npred4<-predict(nline4,newdata=new)
> plot(x,y)
> lines(new$x,npred4,col="gray",lwd=2)


从图中我们可以看出,倒数函数的拟合曲线,确实也效果不太好,而且最后一个点有个想离群点,如果没有最后一个点这个效果应该是更好的。

4.5 模型优化
根据上面的几种模型验证的结果,我手动进行了优化,发现高次函数进行拟合效果会更好。定义公式 y ~ a*x^3+b*x^5+c*x^7+e*d^9+e*x^11 ,其中高次函数定义为11次,并去掉了截距。


> nline5 <- nls(y ~ a*x^3+b*x^5+c*x^7+d*x^9+e*x^11,
+               start = list(a=1,b=2,c=1,d=1,e=1))

# 查看模型结果
> nline5
Nonlinear regression model
  model: y ~ a * x^3 + b * x^5 + c * x^7 + d * x^9 + e * x^11
   data: parent.frame()
         a          b          c          d          e 
 2.959e-01 -2.064e-02  5.837e-04 -7.324e-06  3.367e-08 
 residual sum-of-squares: 2.583

Number of iterations to convergence: 3 
Achieved convergence tolerance: 2.343e-07

# 模型分析
> summary(nline5)

Formula: y ~ a * x^3 + b * x^5 + c * x^7 + d * x^9 + e * x^11

Parameters:
    Estimate Std. Error t value Pr(>|t|)   
a  2.959e-01  5.132e-02   5.767   0.0022 **
b -2.064e-02  6.176e-03  -3.341   0.0205 * 
c  5.837e-04  2.468e-04   2.365   0.0644 . 
d -7.324e-06  3.940e-06  -1.859   0.1222   
e  3.367e-08  2.139e-08   1.574   0.1763   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7187 on 5 degrees of freedom

Number of iterations to convergence: 3 
Achieved convergence tolerance: 2.343e-07

通过查看模型的结果数据,我们可以发现通过T检验相关系数a,b,c是显著,残差的误差平方和(residual sum-of-squares)为2.583,数值是最小的,而且比三次函数增长了2上系数,效果应该算不错的,能够表达拟合的特征。

把观测值与预测值进行可视化展示,以散点表现观测值,以线表示预测值。


> npred5<-predict(nline5,newdata=new)
> plot(x,y)
> lines(new$x,npred5,col="purple3",lwd=2)

从图中我们可以看出,拟合的曲线穿过观测值的区域,并且非常贴近部分数据点,直观上说明效果很好,与数据检验给出的结论是一致的。

4.6 总结

我们把上面的1种线性函数和5种非线性的拟合函数的模型评价指标,整理到一个表格中。其中,三次函数和高次函数的拟合效果是不错的,而其他的几种拟合效果,都是不太好的。判断好坏的标准,就是残差的误差平方和,同时参数的T检验显著,才是比较好的拟合模型。

拟合函数残差平方和参数T检验显著性解释性颜色
一次函数10.95334不显著红色
二次函数9.758不显著蓝色
三次函数3.453显著不好绿色
指数函数8.966不显著不好黄色
倒数函数15.74显著灰色
高次函数2.583显著无法解释紫色

最后,合并所有拟合曲线,画到同一张图上,进行对比。绿色的三次函数拟合曲线,和紫色的高次函数拟合曲线效果最好。

但是,高次函数也带来了的负面问题,一是过拟合,二是业务的可解释性非常差。如果x和y分别带入到业务场景,比如x是货物的重量和y是货物的单价,那么x的11次方,完全是无法理解的,等同于一个黑盒的模型。所以,真正落地项目在使用的时候,要根据业务场景的要求,在准确度,可解释性,复杂度等标准之下,进行权衡。

本文利用最小二乘法,用R语言进行了线性回归和非线性回归的实验。验证了最小二乘法在处理回归问题上是非常有效的,除了回归问题,最小二乘法还能解决分类的问题。最小二乘法作为机器学习的底层算法,是非常值得我们学习和掌握的,明白原理就能掌握一类的高级模型。

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

打赏作者

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

打赏作者