R package 'FNN' (1) カルバックライブラダイバージェンス

R Program Mathematics 
#install.packages("FNN")
library(FNN)

#http://nakhirot.hatenablog.com/entry/20130607/1370602095　
#で理解したKL.divergenceについて、Rの関数(KL.divergence())を分解して理解する

####KL.divergence####
function (X, Y, k = 10, algorithm = c("VR", "brute", "kd_tree","cover_tree"))
#X,Yはmatrixで、XとYの列数は等しくなければならない。
#kはhelpによると、最も近い値を周囲幅いくつまで探すのかを意味する。

{
    #match.argは選択肢を受け取る引数を作るための関数(http://d.hatena.ne.jp/hoxo_m/20111025/p1)
    algorithm <- match.arg(algorithm)

    #inputのマトリックス化
    if (!is.matrix(X))
        X <- matrix(X)
    if (!is.matrix(Y))
        Y <- matrix(Y)

　　#KL divergenceの計算
    n <- nrow(X)
    p <- ncol(X)
    m <- nrow(Y)

    log(m/n) + p * (colMeans(log(knnx.dist(Y, X, k = k, 
    algorithm))) - colMeans(log(knn.dist(X, k = k, 
    algorithm))))

　　 #colMeans(log(knnx.dist())),colMeans(log(knn.dist()))の項はKL-divergence
　　 #(http://nakhirot.hatenablog.com/entry/20130607/1370602095)の
　　 #項-p*logq, p*logpにそれぞれ対応していると想定される。
　　 #引き数X,Yの行数が異なる場合は、その分を調節する。
  　 #Xの行数が大きい場合はその分KL.divergenceの距離を小さくする(log(m/n))
}
###knn.dist###
function (data, k = 10, algorithm = c("cover_tree", "kd_tree",
    "VR", "CR", "brute"))
{
    get.knn(data, k, algorithm)$nn.dist
}
#get.knn関数のnn.dist成分を抽出しているだけ。

##get.knn##
function (data, k = 10, algorithm = c("cover_tree", "kd_tree",
    "VR", "CR", "brute"))
{
    #match.argは選択肢を受け取る引数を作る関数(http://d.hatena.ne.jp/hoxo_m/20111025/p1)
    algorithm <- match.arg(algorithm)

    #データのmatrix化
    if (!is.matrix(data))
        data <- as.matrix(data)
　　 #データが数値(numeric)でない場合はエラーを表示
    if (!is.numeric(data))
        stop("Data non-numeric")

　　 #データにNAが含まれる場合はエラーを表示
    if (any(is.na(data)))
        stop("Data include NAs")


　　 #storage.modeはデータの型を調べる関数(http://ofmind.net/doc/r-intro-lecture)
　　 #classよりは細かい区別(double,integerなど)を返す
    #データ型をintegerからdoubleに変換
    if (storage.mode(data) == "integer")
        storage.mode(data) <- "double"
    
    n <- nrow(data)
    d <- ncol(data)

    #inputの行数が引数k以下である場合は警告を表示
    if (k >= n)
        warning("k should be less than sample size!")

　　 #switchについては、http://nakhirot.hatenablog.com/search?q=switch
    Cname <- switch(algorithm, cover_tree = "get_KNN_cover",
        kd_tree = "get_KNN_kd", VR = "get_KNN_VR", CR = "get_KNN_CR",
        brute = "get_KNN_brute")

    #.CはRからコンパイル済みのCコードを呼び出すコマンド(こちらを参照)
　  #nn.indexはn*k個の数列(integer型), nn.distはn*k個の数列(double型)
    knnres <- .C(Cname, t(data), as.integer(k), d, n, nn.index = integer(n *
        k), nn.dist = double(n * k), DUP = FALSE)

 　 #nn.indexをn行k列の行列に整形
 　 nn.index <- matrix(knnres$nn.index, byrow = T, nrow = n, ncol = k) 
 
　  #nn.distをn行k列の行列に整形
    nn.dist <- matrix(knnres$nn.dist, byrow = T, nrow = n, ncol = k) 
　　
    #引数k(デフォルト10)よりinput行列の行数nが小さい場合はnn.index,nn.listのn～k列をNAとする
    if (k >= n) { nn.index[, n:k] <- NA nn.dist[, n:k] <- NA }
　　
    #nn.index,nn.distを返す
    return(list(nn.index = nn.index, nn.dist = nn.dist)) }

 

###knnx.dist###
function (data, query, k = 10, algorithm = c("cover_tree", "kd_tree",
    "VR", "CR", "brute"))
{
    get.knnx(data, query, k, algorithm)$nn.dist
}
#knn.distと同様、get.knnxのnn.dist成分を返すだけの関数。

##get.knnx (略):get.knnと同様

【まとめ】

  • KL.divergence()の最下層にあるget.knn, get.knnxはR上からCのコードを呼び出して実行しているため、細部まで淀みなく理解するにはC言語の習得が必要。
  • 2つの行列X,YのKL.divergenceの計算時には、XとYの列数が等しくなければならない(任意の2つの観測値の距離を計算するため)。各行は1つの観測に対応するものと認識。
  • get.knn, get.knnxで出力されるnn.indexは各観測(行番号)に"近い"観測(行番号)が1番近いものから順に列番号順に並べられている。
  • nn.distはその場合の各観測間の距離。
  • KL.divergence関数では、get.knn,get.knnxの出力(nn.dist)の対数(log)をとり平均値(colMeans)をとることにより、1番目近い観測とのlog(距離)の平均、2番目に近い観測とのlog(距離)の平均、…を出力している。
  • 最後に、観測数(行数)差異を補う項(log(m/n))が加算されている。
  • 類似のKL.dist関数は、KL.divergence(X,Y)とKL.divergence(Y,X)を足し合わせて算出。これにより、X,Yに関して対称化されるため「距離」と見なすことが出来るようになると認識。

【追記】

KL.divergenceは、1つの目的変数と多数の独立変数の関係が分布で与えられている際に、多数の目的変数と独立変数の観測そ含む多数のサンプル間の関係を調べるために、用いられる。（という理解）