在R中实现二维核密度估计的不同内核

Question

在R中实现二维核密度估计的不同内核

4

我需要帮助理解如何实现一个二维核密度方法，使用各向同性方差和双变量正态核。然而，由于数据位于地球表面上，因此需要使用大圆距离，而不是常规距离。我想在R中复制这个过程，但我无法弄清楚如何使用除简单欧几里得距离之外的距离度量来估计卷积核。是否有人能够编写任意核函数的程序？

- David Manheim

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- David Manheim · Accepted Answer

我最终修改了MASS库中的kde2d函数。如下所示，需要进行一些重大修订。尽管如此，该代码非常灵活，允许使用任意的二维核。在这种情况下，采用rdist.earth()来计算大圆距离，h是所选带宽，以千米为单位，n是要在每个方向上使用的网格点数。rdist.earth需要“fields”库。

虽然可以修改该函数以执行超过2d的计算，但是在更高的维度中，网格会非常快地变大。(现在已经不小了。)

欢迎对优雅性或性能提出评论和建议！

kde2d_mod <- function (data, h, n = 200, lims = c(range(data$lat), range(data$lon))) {
#Data is a matrix: lon,lat for each source. (lon,lat to match rdist.earth format.)
print(Sys.time()) #for timing

nx <- dim(data)[1]
if (dim(data)[2] != 2) 
stop("data vectors have only lat-long data")
if (any(!is.finite(data))) 
stop("missing or infinite values in the data are not allowed")
if (any(!is.finite(lims))) 
stop("only finite values are allowed in 'lims'")
#Grid:
g<-grid(n,lims) #Function to create grid.

#The distance matrix gets large... Can we work around it? YES WE CAN!
sets<-ceiling(dim(g)[1]/10000)
#Allocate our output:
z<-rep(as.double(0),dim(g)[1])

for (i in (1:sets)-1) {
   g_subset=g[(i*10000+1):(min((i+1)*10000,dim(g)[1])),]
   a_matrix<-rdist.earth(g_subset,data,miles=FALSE)

   z[(i*10000+1):(min((i+1)*10000,dim(g)[1]))]<- apply( #Here is my kernel...
    a_matrix,1,FUN=function(X)
    {sum(exp(-X^2/(2*(h^2))))/(2*pi*nx)}
   )
rm(a_matrix)
}

print(Sys.time())
#Un-transpose the final data.
z<-t(matrix(z,n,n))
dim(z)<-c(n^2,1)
z<-as.vector(z)
return(z)
}

这里的关键点是内循环中可以使用任何内核；缺点是这会被评估在网格点上，因此需要高分辨率的网格来运行；FFT会很好，但我没有尝试过。

网格函数：

grid<- function(n,lims) {
num <- rep(n, length.out = 2L)
gx <- seq.int(lims[1L], lims[2L], length.out = num[1L])
gy <- seq.int(lims[3L], lims[4L], length.out = num[2L])

v1=rep(gy,length(gx))
v2=rep(gx,length(gy))
v1<-matrix(v1, nrow=length(gy), ncol=length(gx))
v2<-t(matrix(v2, nrow=length(gx), ncol=length(gy)))
grid_out<-c(unlist(v1),unlist(v2))

grid_out<-aperm(array(grid_out,dim=c(n,n,2)),c(3,2,1) ) #reshape
grid_out<-unlist(as.list(grid_out))
dim(grid_out)<-c(2,n^2)
grid_out<-t(grid_out)
return(grid_out)
}

你可以使用image.plot函数来绘制值，其中v1和v2矩阵为x，y坐标：

kde2d_mod_plot<-function(kde2d_mod_output,n,lims) ){
 num <- rep(n, length.out = 2L)
 gx <- seq.int(lims[1L], lims[2L], length.out = num[1L])
 gy <- seq.int(lims[3L], lims[4L], length.out = num[2L])

 v1=rep(gy,length(gx))
 v2=rep(gx,length(gy))
 v1<-matrix(v1, nrow=length(gy), ncol=length(gx))
 v2<-t(matrix(v2, nrow=length(gx), ncol=length(gy)))

 image.plot(v1,v2,matrix(kde2d_mod_output,n,n))
 map('world', fill = FALSE,add=TRUE)
}