模拟 LGCP kppm 时指定一组点数
Specifying a set number of points when simulating a LGCP kppm
似乎有一些方法可以在模拟某些 ppm 对象时使用一组点数(例如,通过操纵 rmhcontrol 选项),但我没有看到 LGCP 的类似选项。
查看文档和源代码后,看起来 rLGCP 的点数是由 lambda 和正在模拟的 window 区域的集成控制的(根据 [=10 中的文档) =] 和对 mu=integrate(lambda)
的调用)。然后在对 rpois
的调用中将其用作“mu”参数,它确定模拟中使用的点数 (nn=rpois(mu,nsims)
)。 因此,目前可用于直接控制点数的唯一方法似乎是编辑各种函数调用和随后的“nn”变量以采用用户定义的常量值,而不是随机抽取基于拉姆达。对吗?
或者,我正在考虑使用随机细化来解决问题,这样我会先进行模拟,然后细化到所需的点数。
这是一个“条件模拟”的问题——你想生成一个点过程模型的模拟实现,条件是点数正好等于n。
这项任务对于泊松点过程模型来说很容易,对于吉布斯点过程模型来说相对简单,这是由 ppm
个对象表示的两种类型的模型。
但是,对于 Cox 过程(包括对数高斯 Cox 过程 LGCP),条件模拟要困难得多。有解决方案,但它们要么在数学上很复杂,要么计算量很大。所以我们还没有在spatstat
.
中实现它们
(当然,您可以修改代码,使其生成具有指定点数的模式,但这不是 LGCP 模型的有效条件模拟)。
但是如果你真的非常需要做条件模拟,试试下面的暴力破解方法。这里 model
是 class 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 过程的条件模拟(固定点数)现已在最新版本的包 spatstat.core
中的 simulate.kppm
中实现。参数n.cond
指定固定点数
似乎有一些方法可以在模拟某些 ppm 对象时使用一组点数(例如,通过操纵 rmhcontrol 选项),但我没有看到 LGCP 的类似选项。
查看文档和源代码后,看起来 rLGCP 的点数是由 lambda 和正在模拟的 window 区域的集成控制的(根据 [=10 中的文档) =] 和对 mu=integrate(lambda)
的调用)。然后在对 rpois
的调用中将其用作“mu”参数,它确定模拟中使用的点数 (nn=rpois(mu,nsims)
)。 因此,目前可用于直接控制点数的唯一方法似乎是编辑各种函数调用和随后的“nn”变量以采用用户定义的常量值,而不是随机抽取基于拉姆达。对吗?
或者,我正在考虑使用随机细化来解决问题,这样我会先进行模拟,然后细化到所需的点数。
这是一个“条件模拟”的问题——你想生成一个点过程模型的模拟实现,条件是点数正好等于n。
这项任务对于泊松点过程模型来说很容易,对于吉布斯点过程模型来说相对简单,这是由 ppm
个对象表示的两种类型的模型。
但是,对于 Cox 过程(包括对数高斯 Cox 过程 LGCP),条件模拟要困难得多。有解决方案,但它们要么在数学上很复杂,要么计算量很大。所以我们还没有在spatstat
.
(当然,您可以修改代码,使其生成具有指定点数的模式,但这不是 LGCP 模型的有效条件模拟)。
但是如果你真的非常需要做条件模拟,试试下面的暴力破解方法。这里 model
是 class 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 过程的条件模拟(固定点数)现已在最新版本的包 spatstat.core
中的 simulate.kppm
中实现。参数n.cond
指定固定点数