自分のための「RでPAMの応用」

この数日、調査やら、学会やらでブログをUPできなかった。

今回の調査では、極めて現代社会史的な事を調べるために、
私と同じ人物を追いかけている人物に会いに行った。
結果として、非常に面白い話が聞けたので収穫は多かった。

ところで、私が会いに行った方というのが、なんと、
UNIXの話やアセンブリの話もよく知っておられて、
80年代のコンピュータ技術について大いに盛り上がった。

いやいや、改めて、文系と理系の区別がどうでも良いことが解った。
結局のところ「やる」か、「やらない」か、の問題であって、
「文系出身だから…。」というのは、言い訳にしか過ぎないのである。

やはり、文系であっても機械の一つや二つ、作れないといけない。
私は、まだまだ、入門したての素人ではあるが、とりあえず、
成層圏まで上がれる何かの観測機器開発に挑戦してみよう。

何やら、出来る気がしてきた。いや、出来ないかもしれないが。
出来なくてもいいや。これで研究助成金を狙うわけでもない。
ただ、宇宙関連のことは自分でもやってみたい。

学会では、ハンズオン・セッションなるものに参加した。
国内学会でこの手のイベントに参加したのは初めてだったので、
それなりに楽しめたと思う。面白い企画があれば、今後も参加してみよう。

ふむ。何やら、今回は、全くどうでも良い話ばかりしている。
ここからが、本日のブログの本題。

今回は、Rのハンズオン・セッションにも参加したのであるが、
そのチュートリアルでは、さり気なく「流行り」のクラスタ分析を行なっていた。
それは、「Partition Around Medoids（PAM）」というもの。

結局、時間が足りなくなって、PAMの話はスルーであったのだが、
素人向けの講習会で、PAMを使うと言うのは正直驚いた。良い意味で。
もっとも、この方法がベストであるというわけではないが。

詳しい話は、「文系のための多次元データ分析」で書く予定であるが、
PAMやK-Means法などの階層的クラスタリングには、
妥当な分類基数を知る方法が存在している。

その方法には、シルエット幅というのを利用するものがあり、
ある分類の誤判別の程度を使って、分類基数の候補を絞りこむ。
簡便な方法であるが、Rの関数として定義されていない。

ということで、学会からの帰りの新幹線の中でコーディング。

以下の関数を用いると、分類基数を推薦してくれる。
なお、クラスタ分析を複数回繰り返すので、
巨大なデータにはあまりには向かない。

以下、ソースコード。

---

cis.estimateCluster <- function(x, from=2, to=10, ...){
  # Load a required package.
  require(cluster)

  # Get number of clusters for checking.
  n <- to-from + 1

  # Create a vector to store the result of silhouette width.
  x.sil <- vector(length=n+1)
  names(x.sil) <- c(seq(from, to),"")

  # Initialyse a step counter.
  cnt <- 1

  # Get a silhouette width in each cluster.
  for(i in from:to){
    sil <- pam(x, k=i)$silinfo$avg.width
    x.sil[cnt] <- mean(sil)
    cnt <- cnt + 1
  }

  # Get the highest value of the silwidth.
  k <- which(x.sil == max(x.sil))
  res <- as.numeric(names(x.sil)[k])

  # Sort with average silhoutte width.
  x.srt <- sort(x.sil, decreasing = TRUE)

  # Create labels and legends for plot.
  labels <- matrix(rep(c("",0), n),ncol=2, byrow=TRUE)
  legend <- vector(length=5)

  # Suggest the numbers of culsters.
  title <- "Suggested numbers for clustering"
  print(title)
  for (i in 1:5){
    # Get number of clustering and average silhouette width.
    x.cls.num <- as.numeric(names(x.srt)[i]) # The number of sugested clusters.
    x.cls.val <- as.numeric(round(x.srt[i], 3)) # Average silhoutte width.

    # Create a label used for legend.
    legend[i] <- paste(names(x.srt)[i], " clusters with ",
                 as.character(x.cls.val), " average silhoutte width.")

    # Create a label used for plot.
    labels[x.cls.num, 1] <- as.character(x.cls.num)
    labels[x.cls.num, 2] <- x.cls.val

    # Output the result.
    print(legend[i])
  }

  # Plot a graph for average silhouette width in each cluster.
  plot(x.sil, type="l", xaxt="n", xlab="", ylab="", col="red"); par(new=TRUE)
  plot(x.sil, type="h", xaxt="n", xlab="", ylab="", lty=2, col="grey")
  title(xlab="Number of cluster", ylab="Average silhoutte width")
  title(main="Average silhoutte width in each class")
  axis(1, labels=names(x.sil), at=seq(1, n+1))
  text(seq(0, n-1), as.numeric(labels[,2]), labels[,1])
  legend("topright", legend, title=title)

  # Return value.
  return(x.sil)
}