density-estimation.R

source('~nikolas/bin/FCS/hdr.R')

library(ggplot2)
library(mvtnorm)

sigma<- function(v, r, p)
{
 V<- matrix(r^2, ncol=p, nrow=p)
   diag(V)<- 1
 V*v
}

# 2 dimensional
X<- data.frame(rmvnorm(1000, mean=rep(0, 2), sigma(1, .5, 2)))
real.dens <- dmvnorm(X, mean=rep(0, 2), sigma(1, .5, 2))

X <- kde2D(X[,1:2],.1)
X$knn.dens <- knn.dens(X,k=10)

ggplot(X, aes(x=X1,y=X2,z=z.knn)) + geom_point() + stat_contour()

ggplot(X, aes(x=X1,y=X2,z=dens)) + geom_point() + geom_density2d()

plot(real.dens, kde2D(X[,1:2],.5)[,3])
plot(real.dens, knn.dens(X,k=100))


# N dimensional
# k neigbours
benchmark.knn.dens <- function(k,N) {
    X <- rmvnorm(10000, mean=rep(0, N), sigma(1, .5, N))
    real.dens <- dmvnorm(X, mean=rep(0, N), sigma(1, .5, N))
    est.dens <- knn.dens(X,k=10)
    x <- real.dens/sum(real.dens)
    y <- est.dens/sum(est.dens)
    lim <- range(c(x,y))
    #plot(x, y, xlab='real density', ylab='knn density', xlim=lim, ylim=lim)
    #abline(b=1,a=0, lwd=2)
    #abline(line(y,x))
    #return(coef(line(y,x))[[2]])
    return(cor(x,y))
}

r.squared <- sapply( 2:100, function(n) benchmark.knn.dens(k=10, n) )
plot(2:100,r.squared, type='l', xlab='dimensions')

set.seed(1234)
r.squared <- sapply( 2:100, function(k) benchmark.knn.dens(k=k, 2) )
plot(2:100,r.squared, type='l', xlab='# of neighbours')