Реализация другого ядра для двумерной оценки плотности ядра в R

В итоге я модифицировал функцию kde2d из библиотеки MASS. Как показано ниже, потребовалась значительная доработка. Тем не менее, код очень гибкий, что позволяет использовать произвольное двумерное ядро. (rdist.earth () использовался для расстояния по большому кругу, h - это выбранная ширина полосы, в данном случае в км, а 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)
}

Ключевым моментом здесь является то, что во внутреннем цикле можно использовать любое ядро; Обратной стороной является то, что это оценивается в точках сетки, поэтому для этого требуется сетка с высоким разрешением; БПФ было бы здорово, но я не пробовал.

Функция сетки:

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)
}