`density.ppp` 比较 `kernel = "gaussian"` 和 Gaussian 内核重建结果的问题

Question

我目前正在集中使用 density.ppp 函数，用我自己设计的不同内核函数调用它。对于我的项目，我还需要 'rebuild' 已经可用的内核函数，例如 "gaussian" 和 "quartic" kernel.

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

我发现对于高斯内核，这与对其他内置内核做同样的事情有很大不同，例如kernel = "quartic".

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

我有两个主要问题：

如果使用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)。我说得对吗？

不知何故，我暂时确信 kernel = "gaussian" 与 Gaussian-kernel-rebuild 的结果也不同，因为在 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

comment 意味着对于 kernel = "gaussian" 进行了更改以避免这种不必要的计算，对吗？那为什么其他内核还在那里呢？

Answer 1

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

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

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

`density.ppp` 比较 `kernel = "gaussian"` 和 Gaussian 内核重建结果的问题

Issues in `density.ppp` comparing results for `kernel = "gaussian"` and Gaussian kernel rebuild by user

r

spatstat

kernel-density