绝大部分代码来源: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) #输出聚类后类间距离
整理后的输出结果:
从左至右依次为:4列输入的x(并用颜色划分了聚类结果,数据太长了后面没放);输出的类间距离di;输出的分割点对照表j(具体查表方法见下文,虽然可能大概率用不上,但是还是啰嗦一嘴)
补充:
此表为Fisher最优分割法计算出的分割点对照表,需从下向上、从右向左看。举例:82行分为5类时的分割点为76,即第五类为{76至82},随后跳转到75行,75行分5类时的分割点为44,即第五类为{44至75},以此类推完成分割。
本人菜鸡一个,水平有限,有书写错误的地方还望各位指出更正。