# カーネル密度分布推定

• １次元、２次元プロットから密度分布推定をする
• 分布(たとえば正規分布)の重なりとして推定する
```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)
```