Kmeans法のソースを見つつ

Cの勉強もかねつつ。kmeans.cとkmeans.Rから最低限必要そうなところを引っ込抜いてきて、コメントを適度に埋めつつやってみました。

Cのソース。100行くらい。短かい。

#include <R.h>
#include <Rinternals.h>

void mykmeans(double *x, int *pn, int *pp, double *cen, int *pk, int *cl,
		  int *pmaxiter, int *nc, double *wss) 
{
  int n = *pn, k = *pk, p = *pp, maxiter = *pmaxiter;
  int iter, i, j, c, it, inew = 0;
  double best, dd, tmp;
  Rboolean updated;
  for(i = 0; i < n; i++) cl[i] = -1;
  /* cl[i]はi番目のデータが第何クラスタに所属しているかを表わす */
  for(iter = 0; iter < maxiter; iter++) {
    Rprintf("%d回目のiterationです。\n",iter+1);
    updated = FALSE;
    for(i = 0; i < n; i++) {
      /* データに対するループ */
      /* それぞれの点で、最も近いクラスタの中心(重心)を見つける */
      best = R_PosInf;
      for(j = 0; j < k; j++) {
	/* クラスタ数に対するループ */
	dd = 0.0;
	for(c = 0; c < p; c++) {
	  /* 次元に対するループ。点x_iにおける、各次元でのj番目のクラスタの中心との二乗和を取っている */
	  tmp = x[i+n*c] - cen[j+k*c];
	  dd += tmp * tmp;
	}
	if(dd < best) {
	  /* x_iと一番近いクラスタとの距離に更新 */
	  best = dd;
	  inew = j+1;
	  Rprintf(" %3d番目のデータ(",i+1);
	  for(c = 0; c < p; c++) {
	    if(c < p - 1){
	      Rprintf("% .2f, ",x[i+n*c]);
	    }else{
	      Rprintf("% .2f",x[i+n*c]);
	    }
	  }
	  Rprintf(")がクラスタ%dに移動しました。\n",j+1);
	}
      }
      if(cl[i] != inew) {
	updated = TRUE;
	cl[i] = inew;
      }
    }
    if(!updated) break;
    /* 各クラスタの中心の点を更新 */
    for(j = 0; j < k*p; j++) cen[j] = 0.0;
    for(j = 0; j < k; j++) nc[j] = 0;
    for(i = 0; i < n; i++) {
      it = cl[i] - 1; nc[it]++;
      for(c = 0; c < p; c++) cen[it+c*k] += x[i+c*n];
    }
    for(j = 0; j < k*p; j++) cen[j] /= nc[j % k];
  }

  *pmaxiter = iter + 1;
  for(j = 0; j < k; j++) wss[j] = 0.0;
  for(i = 0; i < n; i++) {
    it = cl[i] - 1;
    for(c = 0; c < p; c++) {
      tmp = x[i+n*c] - cen[it+k*c];
      wss[it] += tmp * tmp;
    }
  }
}

Rのコード。

dyn.load("/tmp/hoge.so")

mykmeans <- function(x, centers, iter.max = 10) {
  cluster <- function() {
    .C("mykmeans", as.double(x), as.integer(m),
       as.integer(ncol(x)),
       centers = as.double(centers), as.integer(k),
       c1 = integer(m), iter = as.integer(iter.max),
       nc = integer(k), wss = double(k))
  }
  x <- as.matrix(x)
  m <- nrow(x)
  k <- centers
  centers <- x[sample(1 : m, k), , drop = FALSE]
  Z <- cluster()
  centers <- matrix(Z$centers, k)
  dimnames(centers) <- list(1:k, dimnames(x)[[2]])
  cluster <- Z$c1
  out <- list(cluster = cluster, centers = centers, withinss = Z$wss,
              size = Z$nc)
  class(out) <- "mykmeans"
  out
}

print.mykmeans <- function(x, ...) {
  cat("K-means clustering with ", length(x$size), " clusters of sizes ",
      paste(x$size, collapse=", "), "\n", sep="")
  cat("\nCluster means:\n")
  print(x$centers, ...)
  cat("\nClustering vector:\n")
  print(x$cluster, ...)
  cat("\nWithin cluster sum of squares by cluster:\n")
  print(x$withinss, ...)
  cat("\nAvailable components:\n")
  print(names(x))
  invisible(x)
}

n <- 100
x <- rbind(matrix(rnorm(n, sd = 0.3), ncol = 2),
           matrix(rnorm(n, mean = 1, sd = 0.3), ncol = 2))
colnames(x) <- c("x", "y")
cl <- mykmeans(x, 2)

plot(x, col = cl$cluster)
points(cl$centers, col = 1:5, pch = 8)