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

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

有序样品聚类介绍

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

S-Plus程序初步说明

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

一共两个函数——

dstc(x, i, j)：即distance，计算两个样品 $x_i$ 和 $x_j$ 之间的距离，参数$x$为样品矩阵 $(x_1,x_2,\ldots,x_n)^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

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

谢益辉 2005-10-01

有序样品聚类介绍

S-Plus程序初步说明