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

Question

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

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

Answer 1

我最终修改了MASS库中的kde2d函数。需要进行一些重大修订，如下所示。也就是说，代码非常灵活，可以使用任意2-D内核。 （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绘制值，与V1和你的X V2矩阵，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) 
}