单置换检验，轮廓检验及其他程序

#单置换检验生成矩阵
#!perl
use warnings;
use strict;

my $n = $ARGV[1];

open FA, $ARGV[0] or die $!;
while(<FA>)
{
	chomp;
	my ($id, @tmp) = split;
	my @s = @tmp[0..$n-1];
	my @e = @tmp[$n..$#tmp];
	my $n = 0;
	for(my $i = 0; $i < @s; $i ++)
	{
		my @st = @s;
		my @et = @e;
		for(my $j = 0; $j < @e; $j ++)
		{
			my $s = $st[$i];
			$st[$i] = $et[$j];
			$et[$j] = $s;
			print $id."-$n\t".join("\t", @st)."\t".join("\t", @et)."\n";
			$n ++;
		}
	}
}

#轮廓检验
args<-commandArgs(T)
library(fpc)
mydat = read.table(args[1], header = F, sep="\t", row.names=1)
va <-rep(0,nrow(mydat))
for(i in 1:nrow(mydat))
{
	stats <- cluster.stats(dist(as.numeric(mydat[i,])), c(rep(1,16), rep(2,16)))
	va[i] <- stats$avg.silwidth
}
out<-cbind(rownames(mydat), va)
write.table(out, file=paste(args[1], ".silhouette", sep = ""), quote=FALSE, col.names=F, sep="\t", row.names=F)

# 轮廓系数与kmeans的关系，转载别人的程序
K <- 2:8
round <- 30 # 每次迭代30次，避免局部最优
rst <- sapply(K, function(i){
    mean(sapply(1:round,function(r){
        result <- kmeans(dat, i)
        stats <- cluster.stats(dist(dat), result$cluster)
        stats$avg.silwidth
    }))
})
plot(K,rst,type='l',main='轮廓系数与K的关系', ylab='轮廓系数')

# 评估聚类结果，转载别人的程序
old.par <- par(mfrow = c(1,2))
k = 2 # 根据上面的评估 k=2最优
clu <- kmeans(dat,k)
mds = cmdscale(dist(dat,method="euclidean"))
plot(mds, col=clu$cluster, main='kmeans聚类 k=2', pch = 19)
plot(mds, col=iris$Species, main='原始聚类', pch = 19)
par(old.par)

kmeans最佳实践

1. 随机选取训练数据中的k个点作为起始点

2. 当k值选定后，随机计算n次，取得到最小开销函数值的k作为最终聚类结果，避免随机引起的局部最优解

3. 手肘法选取k值：绘制出k--开销函数闪点图，看到有明显拐点（如下）的地方，设为k值，可以结合轮廓系数。

4. k值有时候需要根据应用场景选取，而不能完全的依据评估参数选取。