# Kraby, kobiety Pima i irysy klasteryzujemy alg. k-means (k=2,...,50)
# oraz alg. hierarch. klasteryzacji aglomeracyjnej. Nastepnie szacujemy 
# "separowalnosc" danych na k-czesci, czyli najprostrza miare istotnosci
# podzialu. Dla porownania z historycznymi danymi grupujemy iid
# z rozkladu norm, unif i bardzo niesymetrycznego rozkladu gamma. 
#
# Separowalnosc dla metody k-means liczy sie nastepujaco. Dla macierzy 
# danych X mamy podzial calkowitego rozrzutu danych: 
# tr(cov(X)) = tr(cov(W)) + tr(cov(B))
# Lewa strona nie zalezy od podzialu ani od liczby klastrow. Alg. k-means
# minimalizuje tr(cov(W)) po podzialach na k czesci, jego wynik oznaczany
# dalej w_k zalezy tylko od k. Powyzsza rownosc mozna zapisac w skrocie:
# t = w_k + b_k. Separowalnosc to z def sep_k = b_k/t = 1 - w_k/t.
#
# Separowalnosc podzialu na k-czesci dla metody hierarch. klasteryzacji 
# aglomeracyjnej liczy sie z wysokosci dendrogramu: 
# sep_k = 1 - h_k/h_1, gdzie h_k to minimalna wysokosc dla podzialu
# na k-czesci.
#
# W obu przypadkach separowalnosc jest rosnaca funkcja k przypominajaca 
# krzywe ROC. Mamy sep(1)=0, sep(n)=1, gdzie n - liczba wierszy X. 
# Podobnie jak dla krzych ROC do wyboru optymalnego k mozemy stosowac 
# "regule kciuka", czyli wybierac takie k, ze (k,sep_k) jest najblizszy (0,1).
# Regula czasami zawodzi...
 
library(MASS); X=list(); n=1000
X[[1]]=matrix(rnorm(7*n),n,7)
X[[2]]=matrix(runif(7*n),n,7)
X[[3]]=matrix(rgamma(6*n,2),n,6)
Xoryg=na.omit(rbind(Pima.te,Pima.tr))
X[[4]]=as.matrix(Xoryg[,-8])
X[[5]]=as.matrix(crabs[,4:8])
X[[6]]=as.matrix(iris[,-5])
sepK=matrix(NA,50,6); sepH=sepK
for(k in 1:6) {
  X0=X[[k]]; n0=nrow(X0)
  t0=sum(diag(var(X0)*(n0-1)/n0))
  for( i in 2:50) t0=c(t0,sum(kmeans(X0,i,iter.max=150)$with/n0))
  sepK[,k]=1-t0/t0[1]
  h=hclust(dist(X0),method="average")$h[(n0-1):(n0-50)]
  sepH[,k]=1-h/h[1]
}

#postscript("ileKlastrow.eps", onefile=FALSE, horizontal=FALSE, pointsize=10)
par(pty="s")
tytul=c("7dim norm","7dim unif","6dim gamma","Pima","crabs","iris") 
matplot(1:50, sepK, col=1:6, xlab="Liczba klastrow", ylab="Separowalnosc",
                    main="k-means",pch="*")
legend(20,.4,tytul,col=1:6,pch="*")

x11(); par(mfrow=c(2,3))
for (k in 1:6) {
  plot(sepK[,k],ylim=c(0,1),main=tytul[k])
  points(1:50,sepH[,k],col=2); abline(h=.8,col=3)
  legend(20,.2,c("kmeans","hclust"),col=1:2,pch=1)
}
par(mfrow=c(1,1))
#graphics.off()