有序样品聚类分析 S-Plus 程序

谢益辉 / 2005-10-01

我写了两整天的 S-Plus 程序，放在这里做个纪念吧。

有序样品聚类介绍

顾名思义，有序样品即为按照一定顺序排列的样品，其次序不可打乱。有序样品聚类方法与一般的聚类方法区别在于各样品的 “地位” 不相同，而在一般的聚类方法中，所有样品都是按照距离原则平等地被聚到某一类中，不管原先样品是什么样的顺序。

S-Plus 程序初步说明

1. 程序思想 / 算法：主要难点 / 重点在于两个样品 xi 和 xj 之间距离 $D (i, j)$ 的计算以及递推寻找类内最小距离 e[P(n,k)]=min{e[P(j−1,k−1)]+D(j,n)},k<=j<=n
2. 变量 / 函数意义解释：

一共两个函数 ——

• dstc(x, i, j)：即 distance，计算两个样品 xi 和 xj 之间的距离，参数 $x$ 为样品矩阵 (x1,x2,…,xn)T，每一行 $x_i$ 为一个样品（p 维向量），列为变量；
• OCluster(x, k)：即 Ordinal Cluster，寻找类内最小距离以及相应的类起始点，x 同为样品矩阵，k 为分类数，输出的结果有三个矩阵：样品两两之间的直径矩阵 D (i, j)、类内距离矩阵 e [P (n,k)] 和分类起始位置矩阵。

涉及到的变量名 ——

• iniD：即 initial distance，初始距离设为零，用来累加计算 xi 与 xj 之间的直径；
• n：样品数目（样本量），即为样本矩阵的行数；
• mtrxD：即 matrix of distance，n 维矩阵，用来存放直径 D (i,j)；
• mtrxE：n-2 维矩阵，存放最小类内距离；
• mtrxIndex：即 matrix of Index，n-2 维矩阵，用来存放分类的起始点（样品号）；
• iniMin：即 initial Minimum，初始最小距离，用在寻找最小类内距离的算法中，如果某种分法使得类内距离小于 iniMin，那么将 iniMin 替换为这个距离，把距离存放在矩阵 mtrxE 中，然后记下这个分法的起始点存放在 mtrxIndex 中，接着用下一种分法的距离和 iniMin 来比较，重复上面的过程。
• splt：即 split，k-1 维向量，记录把 n 个有序样品分成 k 类时的 k-1 个分割点的位置。
• rownum：即 row number，矩阵行号，具体解释在后面。
• result：聚类的最终结果，本程序用 “//” 符号将样品序号 1, 2, …, n 进行分割表示最后的分类结果。
• 其它 l, a, b, p, w 等均为循环变量。

程序代码如下（核心程序在于找最小类内距离以及相应的分割点）：

dstc <- function(x, i, j) {
  iniD <- 0
  if (i < j) {
    for (l in i:j) {
      iniD <- iniD + crossprod(x[l, ] -
              apply(matrix(x[i:j, ], j - i + 1), 2, mean))
    }
  }
  return(round(iniD, 3))
}

OCluster <- function(x, k) {
  n <- nrow(x)
  mtrxD <- matrix(, n - 1, n - 1)
  for (a in 2:n) {
    for (b in 1:(a - 1)) {
      mtrxD[a - 1, b] <- dstc(x, b, a)
    }
  }
  print(mtrxD)
  mtrxE <- matrix(, n - 2, n - 2)
  mtrxIndex <- matrix(, n - 2, n - 2)
  for (a in 3:n) {
    iniMin <- dstc(x, 1, 2 - 1) + dstc(x, 2, a)
    mtrxE[a - 2, 1] <- iniMin
    mtrxIndex[a - 2, 1] <- 2
    for (p in 3:a) {
      if (dstc(x, 1, p - 1) + dstc(x, p, a) < iniMin) {
        iniMin <- dstc(x, 1, p - 1) + dstc(x, p, a)
        mtrxE[a - 2, 1] <- iniMin
        mtrxIndex[a - 2, 1] <- p
      }
    }
  }
  for (a in 4:n) {
    for (b in 3:(a - 1)) {
      iniMin <- dstc(x, b, a)
      mtrxE[a - 2, b - 1] <- iniMin
      mtrxIndex[a - 2, b - 1] <- b
      for (p in (b + 1):a) {
        if (mtrxE[p - 3, b - 2] + dstc(x, p, a) < iniMin) {
          iniMin <- mtrxE[p - 3, b - 2] + dstc(x, p, a)
          mtrxE[a - 2, b - 1] <- iniMin
          mtrxIndex[a - 2, b - 1] <- p
        }
      }
    }
  }
  print(mtrxE)
  print(mtrxIndex)
  splt <- NULL
  rownum <- n
  for (p in k:2) {
    splt <- c(mtrxIndex[rownum - 2, p - 1], splt)
    rownum <- mtrxIndex[rownum - 2, p - 1] - 1
  }
  result <- NULL
  for (p in 1:n) {
    result <- cat(result, p, " ")
    for (w in 1:length(splt)) {
      if (p == splt[w] - 1) {
        result <- cat(result, "// ")
      }
    }
  }
  return(result)
}

> x <- matrix(c(9.3, 1.8, 1.9, 1.7, 1.5, 1.3, 1.4, 2, 1.9, 2.3, 2.1), 11, 1)
> OCluster(x, 5)
numeric matrix: 10 rows, 10 columns. 
        [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8] [,9] [,10] 
 [1,] 28.125    NA    NA    NA    NA    NA    NA    NA   NA    NA
 [2,] 37.007 0.005    NA    NA    NA    NA    NA    NA   NA    NA
 [3,] 42.207 0.020 0.020    NA    NA    NA    NA    NA   NA    NA
 [4,] 45.992 0.087 0.080 0.020    NA    NA    NA    NA   NA    NA
 [5,] 49.128 0.232 0.200 0.080 0.020    NA    NA    NA   NA    NA
 [6,] 51.100 0.280 0.232 0.087 0.020 0.005    NA    NA   NA    NA
 [7,] 51.529 0.417 0.393 0.308 0.290 0.287 0.180    NA   NA    NA
 [8,] 51.980 0.469 0.454 0.393 0.388 0.370 0.207 0.005   NA    NA
 [9,] 52.029 0.802 0.800 0.774 0.773 0.708 0.420 0.087 0.08    NA
[10,] 52.182 0.909 0.909 0.895 0.889 0.793 0.452 0.087 0.08  0.02
numeric matrix: 9 rows, 9 columns. 
      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] 
[1,] 0.005    NA    NA    NA    NA    NA    NA    NA    NA
[2,] 0.020 0.005    NA    NA    NA    NA    NA    NA    NA
[3,] 0.087 0.020 0.005    NA    NA    NA    NA    NA    NA
[4,] 0.232 0.040 0.020 0.005    NA    NA    NA    NA    NA
[5,] 0.280 0.040 0.025 0.010 0.005    NA    NA    NA    NA
[6,] 0.417 0.280 0.040 0.025 0.010 0.005    NA    NA    NA
[7,] 0.469 0.285 0.045 0.030 0.015 0.010 0.005    NA    NA
[8,] 0.802 0.367 0.127 0.045 0.030 0.015 0.010 0.005    NA
[9,] 0.909 0.367 0.127 0.065 0.045 0.030 0.015 0.010 0.005
integer matrix: 9 rows, 9 columns. 
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] 
[1,]    2   NA   NA   NA   NA   NA   NA   NA   NA
[2,]    2    4   NA   NA   NA   NA   NA   NA   NA
[3,]    2    5    5   NA   NA   NA   NA   NA   NA
[4,]    2    5    6    6   NA   NA   NA   NA   NA
[5,]    2    5    5    6    6   NA   NA   NA   NA
[6,]    2    8    8    8    8    8   NA   NA   NA
[7,]    2    8    8    8    8    8    8   NA   NA
[8,]    2    8    8   10   10   10   10   10   NA
[9,]    2    8    8   10   11   11   11   11   11
1  // 2  3  4  // 5  6  7  // 8  9  // 10  11

唉，谁要是看懂了俺的破程序我真要请他 / 她吃饭！编得真不容易，虽说只有七十多行，也是挖空心思……

← 大开杀戒！！！ 8844 域名不出我所料 →