样本进行有序聚类（最优切割）

对样本进行有序聚类（最优切割）
这是做一个项目的时候搜到的代码，我向量化了其中的部分代码，整体结构跟原作的一致，由于出处也不知道哪里，到处都搜得到，所以如有侵权，还望告知。
以下为代码：
ocluster = function(datasam, classnum) {
#有序样本聚类，输入datasam为样本数据阵，每一行为一个样本；
#输入classnum为要分的类数
#返回值result1为分类结果示意图
#各类的起始点存在变量breaks中
#输出三个矩阵 ra_dis:距离矩阵 leastlost:最小损失矩阵 classid:分类标识矩阵
#author:banmudi 2010.11

#样本数
datasam <- as.matrix(datasam)
sam_n = dim(datasam)[1]
#子函数，计算i-j个样本组成的类的半径 #####书上管这个叫直径，不过没事，都是衡量类的大小的
############################################################################# part0
#从part0到part1这段是我写的程序，向量化了作者写的代码以提高效率。是part1到part2的替代品。
radi = function(a) {
i <- min(a) ; j <- max(a)
#提取i-j个样本
temp =as.matrix( datasam[i:j, ])
mu = colMeans(matrix(temp,j-i+1))
vec = apply(matrix(temp,j-i+1), 1, function(x) {
x - mu
})
round(sum(apply(matrix(vec,j-i+1), 2, crossprod)),3) ####这里他写的跟算法有出入，但是结果是一样的，咱们还是严谨一些，按照算法来
}
rd_temp <- as.matrix(expand.grid(1:sam_n , 1:sam_n))
colnames(rd_temp) <- NULL
ra_dis <- matrix(apply(rd_temp , 1 , radi) , nrow = sam_n , ncol = sam_n)
############################################################################# part1
# #从part1到part2这段是作者写的程序，不过在计算距离矩阵是没有使用向量化操作
# #在样本多的时候会可能运行速度慢，所以我向量化了这段代码以提高效率
# radi = function(i, j) {
# #提取i-j个样本
# temp =as.matrix( datasam[i:j, ])
# mu = colMeans(matrix(temp,j-i+1))
# vec = apply(matrix(temp,j-i+1), 1, function(x) {
# x - mu
# })
# round(sum(apply(matrix(vec,j-i+1), 1, crossprod)),3) ###东凡#这里他写的跟算法有出入，但是结果是一样的，咱们还是严谨一些，按照算法来
# }
#
#
# #计算距离矩阵
# ra_dis = matrix(0, sam_n, sam_n)
# rownames(ra_dis) = 1:sam_n
# colnames(ra_dis) = 1:sam_n
# for (i in 1:(sam_n - 1)) {
# for (j in (i + 1):sam_n) {
# ra_dis[i, j] = radi(i, j)
# ra_dis[j, i] = radi(i, j)
# }
# }
########################################################################## part2
#最小损失矩阵，行为样本数，列为分类
#leastlost[i,j]表示把1:i样本分成j类对应的最小损失
leastlost = matrix(, sam_n - 1, sam_n - 1)
rownames(leastlost) = 2:sam_n
colnames(leastlost) = 2:sam_n
diag(leastlost) = 0
#round(leastlost,3);

#记录下对应的分类结点
classid = matrix(, sam_n - 1, sam_n - 1)
rownames(classid) = 2:sam_n
colnames(classid) = 2:sam_n
diag(classid) = 2:sam_n

#分成两类时，填写最小损失阵的第一列
leastlost[as.character(3:sam_n), “2”] = sapply(3:sam_n,
function(xn) {
min(ra_dis[1, 1:(xn - 1)] + ra_dis[2:xn, xn])
})
classid[as.character(3:sam_n), “2”] = sapply(3:sam_n, function(xn) {
which((ra_dis[1, 1:(xn - 1)] + ra_dis[2:xn, xn]) == (min(ra_dis[1,
1:(xn - 1)] + ra_dis[2:xn, xn])))[1] + 1
})
#分成j类时，填写最小损失阵的 第二列到最后一列
for (j in as.character(3:(sam_n - 1))) {
#分成j类
leastlost[as.character((as.integer(j) + 1):sam_n), j] = sapply((as.integer(j) +
1):sam_n, function(xn) {
min(leastlost[as.character(j:xn - 1), as.character(as.integer(j) -
1)] + ra_dis[j:xn, xn])
})

classid[as.character((as.integer(j) + 1):sam_n), j] = sapply((as.integer(j) +                                                                1):sam_n, function(xn) {                                                                  a = which((leastlost[as.character(j:xn - 1), as.character(as.integer(j) -                                                                                                                              1)] + ra_dis[j:xn, xn]) == min(leastlost[as.character(j:xn -                                                                                                                                                                                      1), as.character(as.integer(j) - 1)] + ra_dis[j:xn,                                                                                                                                                                                                                                    xn]))[1] + as.integer(j) - 1                                                                })

}

diag(classid) = 2:sam_n

breaks = rep(0, 1, classnum)
breaks[1] = 1
breaks[classnum] = classid[as.character(sam_n), as.character(classnum)]
flag = classnum - 1
while (flag >= 2) {
breaks[flag] = classid[as.character(breaks[flag + 1] -
1), as.character(flag)]
flag = flag - 1
}

print(“distance matrix:”);#cat(“\n”)
print(ra_dis[2:sam_n,1:(sam_n-1)], na.print = “”); #输出距离矩阵
print(“leastlost matrix:”)
print(leastlost[2:(sam_n-1),1:(sam_n-2)], na.print = “”); #输出最小损失矩阵
print(“classid matrix:”)
print(classid[2:(sam_n-1),1:(sam_n-2)], na.print = “”); #输出分类标识矩阵
cat(“\n”)
plot(leastlost[sam_n - 1,] , type = “b” , main = “最小损失函数随分类数变化的趋势图”,
xaxt = “n” , xlab = “分类数” , ylab = “最小损失函数” , col = “blue”)
axis(1, at = 1:(sam_n-1) , labels =1:(sam_n-1) ,las = 0)
print(“result”)
#画一个简单的分类示意图
result1=NULL
for (p in 1:sam_n) {
result1 <- cat(result1,p, ” “)
for (w in 1:length(breaks)) {
if (p == breaks[w] - 1) {
result1 <- cat(result1, “||”)
}
}
if (p == sam_n)
result1= cat(result1, “\n”)
}
}