R Kernel estimation 2次元カーネル推定

2次元データはこちらで例が載っていたので、ここでは3次元データで2次元カーネル推定を試してみる。MASSパッケージのkde2d関数が、ガウス分布のカーネルを用いて滑らかな密度推定を行ってくれる。

サンプルに用いるデータは、位置情報（lat, long)と各点における犯罪の発生件数（freq）という架空データ。（上記リンクから拝借しました）

library(MASS)

#Data preparation
lat <- c(1,2,2,2,2,3,4,4,5,5,6,6,7,8,8,8,9,11,
         13,15,15,16,16,18,21,28)
long <- c(4,1,2,8,12,7,4,7,3,5,13,14,5,1,3,
          7,15,9,1,11,14,3,4,1,4,11)
freq <- c(11,8,9,4,15,21,11,13,4,7,7,11,17,6,
         17,5,8,1,14,2,5,12,7,6,7,3)

dat <- data.frame(matrix(c(lat,long,freq),ncol=3))
colnames(dat) <- c("lat","long","freq")

datrev <- data.frame(lapply(dat,function(i) rep(i,dat[,"freq"])))

#Kernel estimation
d <- kde2d(datrev$lat,datrev$long,
           c(bandwidth.nrd(datrev$lat),bandwidth.nrd(datrev$long)),n=80)

#Visualization of the kernel estimation
par(mfrow=c(3,1))
plot(datrev$lat,datrev$long,xlab="latitude",ylab="longitude")
image(d,xlab="latitude",ylab="longitude")
contour(d,xlab="latitude",ylab="longitude")

#Detection of points give maximum density
n <- which(d$z == max(d$z),arr.ind=TRUE)
(c(d$x[n[1]],d$y[n[2]]))
#result:lat = 3.392405, long = 7.025316

一番上のグラフは、単にinputデータをプロットしたもの（重なっている点はやや濃くなっている）。中央はカーネル推定の結果（黒い方が犯罪件数が多い）、下段はその等高線。