在 R 中为 2D 核密度估计实现不同的核

Posted 2023-02-14

【中文标题】在 R 中为 2D 核密度估计实现不同的核

【英文标题】：Implementing a different Kernel for 2D Kernel Density Estimation in R

【发布时间】：2013-12-05 05:00:45

【问题描述】：

我正在寻求一些帮助，以了解如何实现具有各向同性方差和二元正态核的二维核密度方法，但不是使用典型的距离，因为数据位于地球表面，我需要使用大圆距离。

我想在 R 中复制这个，但我不知道如何为任何内置估计器使用简单欧几里得距离以外的距离度量，并且因为它使用复杂的卷积方法添加内核。有没有人可以编写任意内核？

【参考方案1】：

我最终修改了 MASS 库中的 kde2d 函数。需要进行一些重大修改，如下所示。也就是说，代码非常灵活，允许使用任意二维内核。 （rdist.earth() 用于大圆距离，h 是选择的带宽，在这种情况下，以 km 为单位，n 是要使用的每个方向上的网格点数。rdist.earth 需要“字段”图书馆）

可以修改该函数以执行超过 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 绘制值，其中 x,y 点使用 v1 和 v2 矩阵：

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)

【讨论】：

在某个时间间隔内（以小时为单位），您可以接受您的答案。 （它似乎不是 kde2d 的直接替代品，因为天真地使用 MASS 中的示例运行它并没有成功。我也收到image(grid) Error in image.default(grid) : increasing 'x' and 'y' values expected 的错误）

这不是替代品； MASS 库假定不相关的 X、Y 内核，这仅在它们处理的非常特定的情况下是正确的。此外， image.plot(output,v1,v2) 对我有用，但仅使用网格函数中的 v1、v2 矩阵；我添加了一个新函数的代码来执行此操作。

仍然得到与with(grid[order(grid$x, grid$y), ], image.plot(x,y,z) ) 相同的错误。我想我的问题是正在绘制哪个对象。抱歉这么密集。

试试新功能。绘制 kde2d_mod 的输出，使用 grid$x,grid$y 作为坐标。