カーネル密度分布推定

• １次元、２次元プロットから密度分布推定をする
  • 分布(たとえば正規分布)の重なりとして推定する

library(KernSmooth)
library(spatgraphs)
# linbin2D関数オブジェクトの呼び出し失敗を回避するために次の処理をする
linbin2D <- get("linbin2D", envir=environment(bkde2D))

# 円周状のプロットを作る
  
Npt<-300
x<-y<-rep(0,Npt)
r<-10
for(i in 1:Npt){
	t<-runif(1)*2*pi/1
	r2<-r*(sin(t*200)+10)+rnorm(1)*0.1
	x[i]<-r2*cos(t)
	y[i]<-r2*sin(t)
}
# spatgraphsのオブジェクトを作る
	pp2d<-list(x=x,y=y,n=length(x),window=list(x=range(c(x,y)),y=range(c(x,y))))
	R<-0.2
	k<-1
   e1<-spatgraph(pp2d,"geometric",par=R)
   e2<-spatgraph(pp2d,"knn",par=k)
   e3<-spatgraph(pp2d,"MST")
      A<-spatcluster(e2)

   #par(mfrow=c(1,3))
   #plot(pp2d,main=paste("Geometric,R =",R))
   #plot(e1,pp2d)
   #plot(pp2d,main=paste("k-nn, k =",k))
   #plot(e2,pp2d)
   #plot(A,pp2d,pch=19)
   plot(pp2d, main="Minimum spanning tree")
   plot(e3,pp2d)


 #f1 <- bkde2D(cbind(x, y), bandwidth=c(width.SJ(x), width.SJ(y)),gridsize=c(101,101))
 f1 <- bkde2D(cbind(x, y) ,bandwidth=c(10,10),gridsize=c(101,101))
  contour(f1$x1, f1$x2, f1$fhat)
persp(f1$fhat)
library(rgl)
xy<-expand.grid(f1$x1,f1$x2)
plot3d(xy[,1],xy[,2],f1$fhat)

# bkde2d の中身を確認するために

x<-cbind(x,y)
gridsize = c(51L, 51L)
truncate = TRUE 
bandwidth=c(10,10)
    
    n <- nrow(x)
    M <- gridsize
    h <- bandwidth
    tau <- 3.4
    if (length(h) == 1L) 
        h <- c(h, h)
    #if (missing(range.x)) {
        range.x <- list(0, 0)
        for (id in (1L:2L)) range.x[[id]] <- c(min(x[, id]) - 
            1.5 * h[id], max(x[, id]) + 1.5 * h[id])
    #}
    a <- c(range.x[[1L]][1L], range.x[[2L]][1L])
    b <- c(range.x[[1L]][2L], range.x[[2L]][2L])
    gpoints1 <- seq(a[1L], b[1L], length = M[1L])
    gpoints2 <- seq(a[2L], b[2L], length = M[2L])
    gcounts <- linbin2D(x, gpoints1, gpoints2)
    L <- numeric(2L)
    kapid <- list(0, 0)
    for (id in 1L:2L) {
        L[id] <- min(floor(tau * h[id] * (M[id] - 1)/(b[id] - 
            a[id])), M[id] - 1L)
        lvecid <- 0:L[id]
        facid <- (b[id] - a[id])/(h[id] * (M[id] - 1L))
        z <- matrix(dnorm(lvecid * facid)/h[id])
        tot <- sum(c(z, rev(z[-1L]))) * facid * h[id]
        kapid[[id]] <- z/tot
    }
    kapp <- kapid[[1L]] %*% (t(kapid[[2L]]))/n
    if (min(L) == 0) 
        warning("Binning grid too coarse for current (small) bandwidth: consider increasing 'gridsize'")
    P <- 2^(ceiling(log(M + L)/log(2)))
    L1 <- L[1L]
    L2 <- L[2L]
    M1 <- M[1L]
    M2 <- M[2L]
    P1 <- P[1L]
    P2 <- P[2L]
    rp <- matrix(0, P1, P2)
    rp[1L:(L1 + 1), 1L:(L2 + 1)] <- kapp
    if (L1) 
        rp[(P1 - L1 + 1):P1, 1L:(L2 + 1)] <- kapp[(L1 + 1):2, 
            1L:(L2 + 1)]
    if (L2) 
        rp[, (P2 - L2 + 1):P2] <- rp[, (L2 + 1):2]
    sp <- matrix(0, P1, P2)
    sp[1L:M1, 1L:M2] <- gcounts
    rp2 <- fft(rp)
    sp2 <- fft(sp)
    rp3 <- Re(fft(rp2 * sp2, inverse = TRUE)/(P1 * P2))[1L:M1, 1L:M2]
    rp4 <- rp3 * matrix(as.numeric(rp3 > 0), nrow(rp3), ncol(rp3))
    list(x1 = gpoints1, x2 = gpoints2, fhat = rp4)