kNN-density estimation

  • クラウド情報からエントロピーを推定する話はこちらに書いた。そこで使っていたのは、kNN(k-nearest neighbor)法
  • その背後には、ある点の密度をその点からk番目に近い点までの範囲に、何個の点があるかで計算しようというもの(参考資料)
  • このkNN密度推定をRでやってみよう
    • ある座標xからk番目に近い点までの距離をr(k,x)とする
    • xを中心とした半径r(k,x)の球は、次元がdの空間であれば、半径r(k,x)のd次元球の体積であるから、それはV(r)= \frac{\pi^{d/2}}{\Gamma(d/2+1)} r^dである
    • この球内に、全部でN個の観察点のうちk個が入っているから
    • P(x) \propto \frac{k}{N} \frac{1}{V(r(k,x))}
  • Rでやってみる
    • x軸、y軸の指定の値についてどんな分布をとるか、と、xy平面にどんな分布が取れるか、をみてやろう


library(FNN)

knn.density <- function(z,X,k=10){
	d <- length(z)
	tmp <- t(t(X)-z)
	D <- sqrt(apply(tmp^2,1,sum))
	D.s <- sort(D)
	(1:k)/length(X[,1])/(pi^(d/2)/gamma(d/2+1)*D.s[1:k]^d)
}

knn.density.multi <- function(Z,X,k=10){
	ret <- matrix(0,length(Z[,1],k))
	for(i in 1:length(Z[,1])){
		ret[i,] <- knn.density(Z[i,],X,k=k)
	}
	ret
}
# m次元データを作る
m <- 2
X <- rbind(cbind(rnorm(200,0,1),rnorm(200,0,5)),cbind(rnorm(300,5,2),rnorm(300,4,1)))
L <- 100
xy <- seq(from=min(X),to=max(X),length=L)

est.D <- matrix(0,length(xy),length(xy))
k <- 10
for(i in 1:L){
	for(j in 1:L){
		est.D[i,j] <- knn.density(c(xy[i],xy[j]),X,k)[k]
	}
}

image(xy,xy,est.D)
points(X,cex=0.1)
contour(xy,xy,est.D,add=TRUE)
par(mfcol=c(1,2))
matplot(est.D,type="l")
matplot(t(est.D),type="l")
par(mfcol=c(1,1))