Saturday, December 13, 2008
A couple of people have asked me to post k-means R code. Here is
Yair's very nice code:
euc <- function(y, x, n, k)
return(sqrt(apply((x - matrix(rep(y,n), nrow=n, ncol=k, byrow=T))^2, 1, sum)))
kmeansyg <- function(x, n_clus, nstart=1)
{
n <- dim(x)[1]
k <- dim(x)[2]
for (i in 1:nstart)
{
### random starting points – make them unique points in the data
repeat
{
starts <- ceiling(runif(n_clus)*n)
if (length(unique(starts)) == n_clus) { break }
}
centers <- x[ceiling(runif(n_clus)*n),]
### initialize clusters so they're unequal, so the while() loop will start
cluster <- rep(0, n)
cluster_old <- rep(1, n)
while (sum(cluster - cluster_old == 0) != n)
{
### save old cluster for while() comparison, and assign new clusters
cluster_old <- cluster
cluster <- apply(apply(centers, 1, euc, x=x, n=n, k=k), 1, which.min)
### compute cluster means
for (j in 1:n_clus)
centers[j,] <- apply(x[cluster == j,], 2, mean)
}
withinss <- 0
for (j in 1:n_clus)
withinss <- (withinss + sum(dist(x[cluster == j,])))/2
withinss <- withinss / length(table(cluster))
test <- list(cluster=cluster, centers=centers, withinss=withinss)
if (i == 1)
answer <- test
if (test$withinss < answer$withinss)
answer <- test
}
return(answer)
}
Yair's very nice code:
euc <- function(y, x, n, k)
return(sqrt(apply((x - matrix(rep(y,n), nrow=n, ncol=k, byrow=T))^2, 1, sum)))
kmeansyg <- function(x, n_clus, nstart=1)
{
n <- dim(x)[1]
k <- dim(x)[2]
for (i in 1:nstart)
{
### random starting points – make them unique points in the data
repeat
{
starts <- ceiling(runif(n_clus)*n)
if (length(unique(starts)) == n_clus) { break }
}
centers <- x[ceiling(runif(n_clus)*n),]
### initialize clusters so they're unequal, so the while() loop will start
cluster <- rep(0, n)
cluster_old <- rep(1, n)
while (sum(cluster - cluster_old == 0) != n)
{
### save old cluster for while() comparison, and assign new clusters
cluster_old <- cluster
cluster <- apply(apply(centers, 1, euc, x=x, n=n, k=k), 1, which.min)
### compute cluster means
for (j in 1:n_clus)
centers[j,] <- apply(x[cluster == j,], 2, mean)
}
withinss <- 0
for (j in 1:n_clus)
withinss <- (withinss + sum(dist(x[cluster == j,])))/2
withinss <- withinss / length(table(cluster))
test <- list(cluster=cluster, centers=centers, withinss=withinss)
if (i == 1)
answer <- test
if (test$withinss < answer$withinss)
answer <- test
}
return(answer)
}