r语言 - "density.ppp"中比较"核 = "gaussian""的结果和用户重建的高斯核的问题

Question

我目前正在深入研究density.ppp函数，用我自己设计的不同内核函数调用它。 对于我的项目，我还需要"重建"已经可用的内核函数，例如"gaussian"和"quartic"内核。

因此，我用kernel = "gaussian"调用density.ppp，并将其结果与使用我自己的内核函数调用density.ppp的结果进行了比较，该函数再现了高斯内核。

我发现对于高斯内核，这与其他内置内核的行为大不相同，例如 kernel = "quartic" .

我知道这是因为内核值的计算方式因所使用的内核而异。对于除"gaussian"之外的所有内核，调用 -one evaluate2Dkernel。另一方面，kernel = "gaussian"的内核值计算被硬编码为 densitypointsEngine 和 second.moment.engine 。

我有两个主要问题：

1. 如果使用at = "points"和leaveoneout = FALSE，

   # ----- contribution from point itself ----------------
   if(!leaveoneout) {
   # add contribution from point itself
   self <- const
   if(!is.null(weights))
     self <- self * weights
   result <- result + self
   

   }

   被执行。

   在此之前，const已设置。对于除 "gaussian" -一，我们有const <- 1/sigma^2.对于kernel = "gaussian"，const是 之后设置为 const <- const/(2*pi) 以便将其用作 硬编码内核中说明函数之前的常量 值计算：

   # constant factor in density computations
   if(is.null(varcov)) {
     const <- 1/sigma^2 
   } else {
     detSigma <- det(varcov)
     Sinv <- solve(varcov)
     const <- 1/sqrt(detSigma)
   }
   if(isgauss) {
     # absorb leading constant in Gaussian density
     const <- const/(2 * pi)
   }
   

   我认为这是不正确的，因为它应该始终为未遗漏的点添加1/sigma^2，但如果kernel = "gaussian"，则会增加1/2*pi*(sigma^2).我说的对吗？

2. 不知何故，我暂时确信我的kernel = "gaussian"与高斯核重建的结果也不同，因为second.moment.engine

   发现

   # set up kernel
   xcol.ker <- xstep * c(0:(nc-1),-(nc:1))
   yrow.ker <- ystep * c(0:(nr-1),-(nr:1))
   kerpixarea <- xstep * ystep
   if(identical(kernel, "gaussian")) {
     if(!is.null(sigma)) {
       densX.ker <- dnorm(xcol.ker, sd=sigma)
       densY.ker <- dnorm(yrow.ker, sd=sigma)
       #' WAS:  Kern <- outer(densY.ker, densX.ker, "*") * kerpixarea
       Kern <- outer(densY.ker, densX.ker, "*")
       Kern <- Kern/sum(Kern)
     } else if(!is.null(varcov)) {
       ## anisotropic kernel
       detSigma <- det(varcov)
       Sinv <- solve(varcov)
       halfSinv <- Sinv/2
       constker <- kerpixarea/(2 * pi * sqrt(detSigma))
       xsq <- matrix((xcol.ker^2)[col(Ypad)], ncol=2*nc, nrow=2*nr)
       ysq <- matrix((yrow.ker^2)[row(Ypad)], ncol=2*nc, nrow=2*nr)
       xy <- outer(yrow.ker, xcol.ker, "*")
       Kern <- constker * exp(-(xsq * halfSinv[1,1]
                                + xy * (halfSinv[1,2]+halfSinv[2,1])
                                + ysq * halfSinv[2,2]))
       Kern <- Kern/sum(Kern)
     } else 
       stop("Must specify either sigma or varcov")
   } else {
     ## non-Gaussian kernel
     ## evaluate kernel at array of points
     xker <- as.vector(xcol.ker[col(Ypad)])
     yker <- as.vector(yrow.ker[row(Ypad)])
     Kern <- evaluate2Dkernel(kernel, xker, yker,
                              sigma=sigma, varcov=varcov, ...) * kerpixarea
     Kern <- matrix(Kern, ncol=2*nc, nrow=2*nr)
     Kern <- Kern/sum(Kern)
   }
   

   仅当内核未"gaussian"时，才使用 kerpixarea 进行乘法。但是经过一些测试，我发现由于Kern <- Kern / sum(Kern)行，是否完成此乘法并不重要。

   但是我想知道为什么它还在那里？这

   #' WAS: Kern <- outer(densY.ker, densX.ker, "*") * kerpixarea

   注释意味着对于kernel = "gaussian"进行了更改以避免这种不必要的计算，对吧？那么，为什么其他内核仍然存在呢？

Answer 1

您在density.ppp的代码中发现了一个错误，其中at="points"、leaveoneout=FALSE 和 kernel 不是默认值。在这种情况下，无法正确计算距离为零的内核值。

该错误将很快在spatstat的开发版本中修复。

关于您的其他问题，恐怕我没有时间解释我们正在用这段代码做什么。它仍在开发中(正如您在评论中看到的那样"Experimental code"(，以使其运行得更快。遗留代码在那里，因此我们可以检查它是否仍然有效。