Fisher有序样本聚类（R语言实现）

缉熙

已于 2022-10-11 23:35:43 修改

阅读量2.8k

点赞数 6

文章标签： r语言聚类机器学习

于 2022-10-11 23:32:47 首次发布

本文链接：https://blog.csdn.net/qq_43513262/article/details/127273247

版权

绝大部分代码来源：https://gutun.plus/study/fisher%E6%9C%89%E5%BA%8F%E6%A0%B7%E6%9C%AC%E8%81%9A%E7%B1%BB%E7%9A%84r%E5%87%BD%E6%95%B0/
作者：GUTUN （侵删）

关于Fisher有序样本聚类

本文侧重点

方法的具体原理这里就不赘述了，实际上，我在应用中发现很多原理的讲解，但是不太好找到代码实现，所以本文的重点在于说明如何在R中利用Fisher有序样本聚类解决实际问题，需要深刻理解原理的请移步其他帖子或书籍。

一个说明代码的例子

在多元统计这本书上有详细的推导证明，但是个人认为若只想快速明白如何实现，不如直接看上面的那个例题，我这里直接贴过来了。这个例题是一维数据，很容易理解，而且其中的符号和代码中的基本一致，供大家对照学习。

在这里插入图片描述

R代码实现

聚类函数：

OrdinalCluster <- function(x, K = 0, D = NULL, pic = F){
  if(!is.data.frame(x)){
    stop("x must be a data frame.")
  }
  n <- nrow(x)
  if(K == 0){
    K <- n
  }
  else{
    if(as.integer(K) != as.numeric(K)) stop("K must be a integer.")
    if(K < 0) stop("K must be positive.")
    if(K > n) stop("K must be less than the sample size.")
  }
  if(is.null(D)){
    D <- function(x, i ,j){
      xp <- as.data.frame(x[i:j,])
      x_mean <- colMeans(xp)
      x_mean_m <- matrix(rep(x_mean, nrow(xp)), nrow(xp), byrow = TRUE)
      return(sum((xp-x_mean_m)^2))
    }
  }
  D_matrix <- matrix(NA, n, n)
  for (i in 1:n) {
    for (j in 1:i) {
      D_matrix[i,j] <- D(x, i, j)
    }
  }
  D_matrix[upper.tri(D_matrix)] <- t(D_matrix)[upper.tri(D_matrix)]
  e_matrix <- matrix(NA, n, K)
  j <- matrix(NA, n, K)
  for (p in 1:K) {
    for (i in p:n) {
      if(p == 1){
        e_matrix[i,p] <- D_matrix[1,i]
        j[i,p] <- 1
      }
      else{
        e_p <- double(i-p+1)
        for (t in p:i) {
          e_p[t-p+1] <- e_matrix[t-1, p-1] + D_matrix[t, i]
        }
        j[i,p] <- which.min(e_p) + p - 1
        e_matrix[i,p] <- e_p[j[i,p]-p+1]
      }
    }
  }
  e_min <- e_matrix[n, ]
  P <- matrix(NA, K, n)
  for (p in 1:K) {
    for (i in p:1) {
      if(i == p){
        P[p,i] <- j[n,p]
      }
      else{
        P[p,i] <- j[P[p,i+1]-1, i]
      }
    }
  }
  result <- list()
  result$D_matrix <- D_matrix
  if(K == 0){
    result$e_min <- e_min
    result$P <- P
  }
  else{
    result$P <- P[K,1:K]
    result$e_matrix <- e_matrix
    result$j<-j
  }
  if(pic == TRUE){
    e_pic <- plot(2:length(e_min),e_min[-1], type = "l", xlab = "K", ylab = "e_min")
    result$e_pic <- e_pic
  }
  return(result)
}

原文中的参数说明
除此之外，因为实际需要我还修改了一下结果输出，可以通过 out$j 输出类似上文例题中表10.5.2括号内的分割点对照表（可能没啥用），可以用来如例题般手动查找，说明分割依据。

附加一个输出函数，还是因为实际需要……我增加了一个输出最终聚类后结果的类间距离的函数，这里用的是使用欧式距离以及最短距离法计算，需要使用其他距离的可以简单更改。

di <- function(out){
  p <- out$P
  d <- out$D_matrix
  di <- list()
  for (i in 1:(length(p)-2)){
    di <- append(di,min(d[p[i]:(p[i+1]-1),p[i+1]:(p[i+2]-1)]))
  }
  di <-append(di,min(d[p[length(p)-1]:(p[length(p)]-1),p[length(p)]:nrow(d)]))
  return(di)
}

应用示例

data <- read.csv('***.csv', header = T) #导入数据
out <- OrdinalCluster(data, K = 5) #使用上文函数
out$D_matrix #输出分类直径
out$P #输出分割点
out$j #输出分割表
di(out) #输出聚类后类间距离