在模拟 LGCP kppm 时指定一组点数

问题描述

在模拟一些 ppm 对象时,似乎有一些方法可以使用一定数量的点(例如,通过操纵 rmhcontrol 选项),但我没有看到 LGCP 的类似选项。

在查看文档和源代码后,看起来 rLGCP 的点数是由 lambda 和正在模拟的窗口区域的集成控制的(每个文档来自 rpoispp调用mu=integrate(lambda))。然后在调用 rpois 时将其用作“mu”参数,该参数确定模拟中使用的点数 (nn=rpois(mu,nsims))。 因此,目前可用于直接控制点数的唯一方法似乎是编辑各种函数调用和随后的“nn”变量以采用用户定义的常量值,而不是基于随机抽取拉姆达。对吗?

或者,我正在考虑使用随机细化的方法,这样我会先模拟,然后细化到所需的点数。

解决方法

这是一个“条件模拟”的问题——你想生成一个点过程模型的模拟实现,条件是点的数量正好等于 n。

这项任务对于 Poisson 点过程模型很容易,而对于 Gibbs 点过程模型则相对简单,这两种模型由 ppm 对象表示。

但是,条件模拟对于 Cox 过程(包括对数高斯 Cox 过程 LGCP)要困难得多。有解决方案,但它们要么在数学上很复杂,要么在计算上很密集。所以我们没有在 spatstat 中实现它们。

(当然,您可以修改代码以使其生成具有指定点数的模式,但这不是 LGCP 模型的有效条件模拟)。

但如果您真的非常需要进行条件模拟,请尝试以下蛮力方法。这里 model 是类 kppm 的拟合模型,nfix 是模式中所需的点数。

bruteBayes <- function(model,nfix,...) {
   stopifnot(is.kppm(model))
   ## stability factor: probability of exactly nfix points for Poisson
   logpn0 <- dpois(nfix,lambda=integrate(intensity(model)),log=TRUE)
   ## generate a large number of realisations of the random driving intensity
   Xlist <- simulate(model,nsim=1000,saveLambda=TRUE,...)
   lamlist <- lapply(Xlist,getElement,name="Lambda")
   ## compute probability of exactly nfix points from each intensity
   mu <- sapply(lamlist,integrate)
   logpn <- dpois(nfix,lambda=mu,log=TRUE)
   ## stabilise
   pna <- exp(logpn - logpn0)
   ## choose one intensity at random,probability proportional to pn
   i <- sample(length(pna),1,prob=pn)
   lami <- lamlist[[i]]
   ## generate 
   Y <- rpoint(nfix,lami)
   return(Y)
}

由于数值下溢,这可能不适用于大型问题。

,

更新:

Cox 过程的条件模拟(保持点数固定)现在在包 simulate.kppm 的最新版本中的 spatstat.core 中实现。参数 n.cond 指定固定点数。