※ 2016.07.27追記
もう少し複雑なサンプルを使った解析と、
ここで使われているダーバージェンスという情報量に関する記事を追加しました。

0. 前振り

言葉の解析をしたい。
ある文章と、ある文章は近いか？という距離を測りたい。

対象となるデータは離散化された「文字」が「順序」を持って並んでいるとする。

ここでは、「遷移確率」を考える。
このページの例をパクると、

いちごじゃむを食べる。
いちごが好き。
いちごのジャムを作る。

この3文の中で、「い」と「ち」と「ご」の出現数は同一であるが、
「い→ち」遷移確率：100%
「ち→ご」遷移確率：100%
「ご→を」「ご→が」「ご→の」遷移確率：33%
という事になっている。

つまり、文節や単語の切れ目では遷移確率が下がる。

文字列から遷移確率を出し、文字列間の距離を定量してクラスタ化できるだろうか？

1. サンプルデータ

5つの文字(A, B, C, D, E)からなる文字列を考える。
文字には色を与えておく。

200文字を横につなげるとこんな感じになるはず。

これを1カラムとして、40サンプルを横に並べる。
ただし、最初の20サンプルと、後の20サンプルでは、語の出現頻度を偏らせておく。

set.seed(1)

dat <- NULL
for(i in 1:20){
  dat <- cbind(dat, sample(c(1:5, 2, 4, 4, 4, 5, 5), 200, replace = T))
}
for(i in 1:20){
  dat <- cbind(dat, sample(c(1:5, 1:5, 2, 3), 200, replace = T))
}

2. 遷移確率マトリックス

横軸はサンプル番号。
縦軸は、「A→A」「A→B」「A→C」..「E→E」という組み合わせ。25種類。
遷移確率=最小値だと白、=最大値だと黒。
この場合、

range(Result)
[1] 0.00 0.75

最初の20サンプルは"D"と"E"が多いので、それに対する遷移に偏っている。
後の20サンプルは、緩やかに"B"と"C"が多くなっているので、パタンが分かれている。

いや、まぁ、どうやっても書けるんですが。
200文字列の、「1~199番目」をdat1に、「2~200番目」をdat2に入れて、

table(dat1, dat2)

とすれば「遷移数」は出るかな、と考える。
ただし、「5種類」が全て出現するとは限らないので、factor(hoge, levels=c(1:5))としているのはそのため。

割り算は、まぁなんとなくapply。
もっと格好良く書けるはず。

Result <- NULL
for(i in 1:ncol(dat)){
dati <- factor(dat[,i], levels=c(1:5))

dat1 <- dati[1:(length(dati)-1)]
dat2 <- dati[2:length(dati)]

a <- function(a){a/table(dati)}
tpi <- apply(table(dat1, dat2), 2, a)

Result <- cbind(Result, c(t(tpi)))
}

3. Jensen-Shannon Divergence

元のデータだと「文字が何番目に現れるか」という情報を持つ。
サンプル1とサンプル2の「N番目」が同じ意味を持つかどうかがよく分からないため、比較していいのかよく分からない。

遷移確率マトリックスに直すと「A→Bの確率」というサンプルを超えて比較可能な情報に変換されているため、比較しやすい。

あとはこの確率のセットをどうやって比較するか、という事になる。
クラスタ化を遠目に見ているので、「対称」で「距離」として使える情報がいい。

エントロピーを規格化して対称になるように補正したJS divergenceを使おう。
これを求めるには、非対称なKullback-Leiber dicergenceを求めて、補正する。
KL divergenceはentropy::KL.pluginで求められる。

もう、for文でいいですよね？

library("entropy")
matJSD <- matrix(0, ncol(Result), ncol(Result))
for(i in 1:ncol(Result)){
  for(k in 1:ncol(Result)){
    a <- Result[,i]
    b <- Result[,k]
    JSD_ik <- try(sum(KL.plugin(b,(a+b)/2), 
                      KL.plugin(a,(a+b)/2))/2, silent=F)
    if(class(JSD_ik)!="try-error"){
      matJSD[i, k] <- JSD_ik
    }
  }
}

警告は、確率０のlogを取ると計算不能になるのを回避するため。

っという事で、

白ければ白いほど、JSDが小さい。

4 clustering

これはもう説明不要でしょう。

plot(hclust(as.dist(matJSD)), hang=-1)

っという事で、メデタシ。