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

谢益辉 2005-10-01

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

有序样品聚类介绍

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

S-Plus程序初步说明

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

一共两个函数——

涉及到的变量名——

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

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

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