Tweaked conditional simulation for Neyman-Scott

Author: baddstats

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: spatstat.model
 Version: 3.3-5.002
-Date: 2025-04-18
+Date: 2025-04-19
 Title: Parametric Statistical Modelling and Inference for the 'spatstat' Family
 Authors@R: c(person("Adrian", "Baddeley", 
                     role = c("aut", "cre", "cph"),
@@ -55,7 +55,7 @@ Authors@R: c(person("Adrian", "Baddeley",
 	     person("Hangsheng", "Wang",
 	            role = "ctb"))
 Maintainer: Adrian Baddeley <[email protected]>
-Depends: R (>= 3.5.0),  spatstat.data (>= 3.1-4), spatstat.univar (>= 3.1-1), spatstat.geom (>= 3.3-4), spatstat.random (>= 3.3-2), spatstat.explore (>= 3.3-0), stats, graphics, grDevices, utils, methods, nlme, rpart
+Depends: R (>= 3.5.0),  spatstat.data (>= 3.1-4), spatstat.univar (>= 3.1-1), spatstat.geom (>= 3.3-4), spatstat.random (>= 3.3-3.005), spatstat.explore (>= 3.3-0), stats, graphics, grDevices, utils, methods, nlme, rpart
 Imports: spatstat.utils (>= 3.1-2), spatstat.sparse (>= 3.1-0), mgcv, Matrix, abind, tensor, goftest (>= 1.2-2)
 Suggests: sm, gsl, locfit, spatial, fftwtools (>= 0.9-8), nleqslv, glmnet, spatstat.linnet (>= 3.2-2), spatstat (>= 3.3)
 Description: Functionality for parametric statistical modelling and inference for spatial data,

NAMESPACE

@@ -110,7 +110,6 @@ export("coef.summary.slrm")
 export("coef.vblogit") 
 export("compareFit")
 export("Concom")
-export("condSimCox")
 export("contour.leverage.ppm")
 export("contour.objsurf")
 export("crosspaircounts")
@@ -305,6 +304,7 @@ export("is.multitype.ppm")
 export("is.multitype.rppm")    
 export("is.multitype.slrm")    
 export("is.poissonclusterprocess")
+export("is.poissonclusterprocess.clusterprocess")
 export("is.poissonclusterprocess.default")
 export("is.poissonclusterprocess.kppm")
 export("is.poissonclusterprocess.zclustermodel")
@@ -326,6 +326,7 @@ export("is.stationary.zgibbsmodel")
 export("Kcom")
 export("killinteraction")
 export("Kmodel")
+export("Kmodel.clusterprocess")
 export("Kmodel.detpointprocfamily")
 export("Kmodel.dppm")
 export("Kmodel.kppm")
@@ -457,6 +458,7 @@ export("parameters.slrm")
 export("parres")
 export("partialModelMatrix")
 export("pcfmodel")
+export("pcfmodel.clusterprocess")
 export("pcfmodel.detpointprocfamily")
 export("pcfmodel.dppm")
 export("pcfmodel.kppm")
@@ -903,6 +905,7 @@ S3method("is.multitype", "msr")
   S3method("is.multitype", "ppm")
 S3method("is.multitype", "rppm")
 S3method("is.multitype", "slrm")
+S3method("is.poissonclusterprocess", "clusterprocess")
 S3method("is.poissonclusterprocess", "default")
 S3method("is.poissonclusterprocess", "kppm")
 S3method("is.poissonclusterprocess", "zclustermodel")
@@ -919,6 +922,7 @@ S3method("is.stationary", "kppm")
 S3method("is.stationary", "ppm")
 S3method("is.stationary", "slrm")
 S3method("is.stationary", "zgibbsmodel")
+S3method("Kmodel", "clusterprocess")
    S3method("Kmodel", "detpointprocfamily")
    S3method("Kmodel", "dppm")
    S3method("Kmodel", "kppm")
@@ -978,6 +982,7 @@ S3method("parameters", "kppm")
 S3method("parameters", "ppm")
 S3method("parameters", "profilepl")
 S3method("parameters", "slrm")
+S3method("pcfmodel", "clusterprocess")
    S3method("pcfmodel", "detpointprocfamily")
    S3method("pcfmodel", "dppm")
    S3method("pcfmodel", "kppm")

R/clusterprocess.R

@@ -0,0 +1,31 @@
+#'
+#'   clusterprocess.R
+#'
+#'   $Revision: 1.1 $ $Date: 2025/04/19 05:17:27 $
+#'
+
+#' clusterprocess() is defined in spatstat.random
+
+#' The following methods are for generics defined in spatstat.model
+
+is.poissonclusterprocess.clusterprocess <- function(model) { TRUE }
+
+pcfmodel.clusterprocess <- function(model, ...) {
+  p <- model$rules$pcf
+  mpar <- model$par.idio
+  other <- model$other
+  f <- function(r) {
+    as.numeric(do.call(p, c(list(par=mpar, rvals=r), other)))
+  }
+  return(f)
+}
+
+Kmodel.clusterprocess <- function(model, ...) {
+  K <- model$rules$K
+  mpar <- model$par.idio
+  other <- model$other
+  f <- function(r) {
+    as.numeric(do.call(K, c(list(par=mpar, rvals=r), other)))
+  }
+  return(f)
+}

R/simulatekppm.R

@@ -3,7 +3,7 @@
 #'
 #'    simulate.kppm
 #'
-#'    $Revision: 1.12 $ $Date: 2023/10/20 11:04:52 $
+#'    $Revision: 1.13 $ $Date: 2025/04/19 05:23:43 $
 
 simulate.kppm <- function(object, nsim=1, seed=NULL, ...,
                           window=NULL, covariates=NULL,
@@ -41,7 +41,7 @@ simulate.kppm <- function(object, nsim=1, seed=NULL, ...,
   if(!is.null(n.cond)) {
     ## fixed number of points
     out <- condSimCox(object, nsim=nsim, seed=NULL, ..., 
-                      window=win, covariates=covariates, 
+                      win=win, covariates=covariates, 
                       n.cond=n.cond, w.cond=w.cond,
                       verbose=verbose, retry=retry, drop=drop)
     out <- timed(out, starttime=starttime)
@@ -144,98 +144,3 @@ simulate.kppm <- function(object, nsim=1, seed=NULL, ...,
   return(out)
 }
 
-condSimCox <- function(object, nsim=1,
-                       ..., window=NULL,
-                       n.cond=NULL, w.cond=NULL,
-                       giveup=1000, maxchunk=100,
-                       saveLambda=FALSE,
-                       verbose=TRUE, drop=FALSE) {
-  if(!inherits(object, c("kppm", "zclustermodel")))
-    stop("object should belong to class 'kppm' or 'zclustermodel'", call=FALSE)
-  shortcut <- isFALSE(object$isPCP)
-
-  w.sim <- as.owin(window)
-  fullwindow <- is.null(w.cond)
-  if(fullwindow) {
-    w.cond <- w.sim
-    w.free <- NULL
-  } else {
-    stopifnot(is.owin(w.cond))
-    w.free <- setminus.owin(w.sim, w.cond)
-  }
-  
-  nremaining <- nsim
-  ntried <- 0
-  accept <- FALSE
-  nchunk <- 1
-  phistory <- mhistory <- numeric(0)
-  results <- list()
-  while(nremaining > 0) {
-    ## increase chunk length
-    nchunk <- min(maxchunk, giveup - ntried, 2 * nchunk)
-    ## bite off next chunk of simulations
-    if(shortcut) {
-      lamlist <- simulate(object, nsim=nchunk,
-                          Lambdaonly=TRUE,
-                          ..., drop=FALSE, verbose=FALSE)
-    } else {
-      Xlist <- simulate(object, nsim=nchunk,
-                        saveLambda=TRUE,
-                        ..., drop=FALSE, verbose=FALSE)
-      lamlist <- lapply(unname(Xlist), attr, which="Lambda", exact=TRUE)
-    }
-    ## compute acceptance probabilities
-    lamlist <- lapply(lamlist, "[", i=w.sim, drop=FALSE, tight=TRUE)
-    if(fullwindow) {
-      mu <- sapply(lamlist, integral)
-    } else {
-      mu <- sapply(lamlist, integral, domain=w.cond)
-    }
-    p <- exp(n.cond * log(mu/n.cond) + n.cond - mu)
-    phistory <- c(phistory, p)
-    mhistory <- c(mhistory, mu)
-    ## accept/reject
-    accept <- (runif(length(p)) < p)
-    if(any(accept)) {
-      jaccept <- which(accept)
-      if(length(jaccept) > nremaining)
-        jaccept <- jaccept[seq_len(nremaining)]
-      naccepted <- length(jaccept)
-      if(verbose)
-        splat("Accepted the",
-              commasep(ordinal(ntried + jaccept)),
-              ngettext(naccepted, "proposal", "proposals"))
-      nremaining <- nremaining - naccepted
-      for(j in jaccept) {
-        lamj <- lamlist[[j]]
-        if(min(lamj) < 0)
-          lamj <- eval.im(pmax(lamj, 0))
-        if(fullwindow) {
-          Y <- rpoint(n.cond, lamj, win=w.sim, forcewin=TRUE)
-        } else {
-          lamj.cond <- lamj[w.cond, drop=FALSE, tight=TRUE]
-          lamj.free <- lamj[w.free, drop=FALSE, tight=TRUE]
-          Ycond <- rpoint(n.cond, lamj.cond, win=w.cond)
-          Yfree <- rpoispp(lamj.free)
-          Y <- superimpose(Ycond, Yfree, W=w.sim)
-        }
-        if(saveLambda) attr(Y, "Lambda") <- lamj
-        results <- append(results, list(Y))
-      }
-    }
-    ntried <- ntried + nchunk
-    if(ntried >= giveup && nremaining > 0) {
-      message(paste("Gave up after", ntried,
-                    "proposals with", nsim - nremaining, "accepted"))
-      message(paste("Mean acceptance probability =",
-                    signif(mean(phistory), 3)))
-      break
-    }
-  }
-  nresults <- length(results)
-  results <- simulationresult(results, nresults, drop)
-  attr(results, "history") <- data.frame(mu=mhistory, p=phistory)
-  if(verbose && nresults == nsim)
-    splat("Mean acceptance probability", signif(mean(phistory), 3))
-  return(results)
-}

R/zclustermodel.R

@@ -11,6 +11,7 @@ zclustermodel <- function(name="Thomas", ..., mu, kappa, scale) {
   if(missing(mu)) stop("The mean cluster size mu must be given")
   if(missing(scale)) stop("The cluster scale must be given")
   rules <- spatstatClusterModelInfo(name)
+  if(!(rules$isPCP)) stop("zclustermodel only supports Neyman-Scott processes")
   par.std <- c(kappa=kappa, scale=scale)
   par.std <- rules$checkpar(par.std, native=FALSE)
   par.idio <- rules$checkpar(par.std, native=TRUE)
@@ -62,26 +63,6 @@ print.zclustermodel <- local({
 
 
 
-pcfmodel.zclustermodel <- function(model, ...) {
-  p <- model$rules$pcf
-  mpar <- model$par.idio
-  other <- model$other
-  f <- function(r) {
-    as.numeric(do.call(p, c(list(par=mpar, rvals=r), other)))
-  }
-  return(f)
-}
-
-Kmodel.zclustermodel <- function(model, ...) {
-  K <- model$rules$K
-  mpar <- model$par.idio
-  other <- model$other
-  f <- function(r) {
-    as.numeric(do.call(K, c(list(par=mpar, rvals=r), other)))
-  }
-  return(f)
-}
-
 intensity.zclustermodel <- function(X, ...) {
   X$par.std[["kappa"]] * X$mu
 }
@@ -121,8 +102,6 @@ reach.zclustermodel <- function(x, ..., epsilon) {
   2 * clusterradius(x, ..., thresh=thresh)
 }
 
-is.poissonclusterprocess.zclustermodel <- function(model) { TRUE }
-
 simulate.zclustermodel <- function(object, nsim=1, ..., win=unit.square()) {
   with(object, {
     switch(name,
@@ -161,18 +140,33 @@ simulate.zclustermodel <- function(object, nsim=1, ..., win=unit.square()) {
                           ...),
                      clustargs))
          },
-         LGCP = {
-           do.call(rLGCP,
-                   resolve.defaults(
-                     list(kappa=par.std[["kappa"]],
-                          scale=par.std[["scale"]],
-                          mu=mu,
-                          win=win,
-                          nsim=nsim,
-                          ...),
-                     clustargs))
-         },
-         stop(paste("Unrecognised model name", sQuote(object$name)),
+         stop(paste("Unrecognised cluster process model name",
+                    sQuote(object$name)),
               call.=FALSE)
          )})
 }
+
+#' The following methods are for generics defined in spatstat.model
+
+is.poissonclusterprocess.zclustermodel <- function(model) { TRUE }
+
+pcfmodel.zclustermodel <- function(model, ...) {
+  p <- model$rules$pcf
+  mpar <- model$par.idio
+  other <- model$other
+  f <- function(r) {
+    as.numeric(do.call(p, c(list(par=mpar, rvals=r), other)))
+  }
+  return(f)
+}
+
+Kmodel.zclustermodel <- function(model, ...) {
+  K <- model$rules$K
+  mpar <- model$par.idio
+  other <- model$other
+  f <- function(r) {
+    as.numeric(do.call(K, c(list(par=mpar, rvals=r), other)))
+  }
+  return(f)
+}
+

inst/doc/packagesizes.txt

@@ -19,4 +19,4 @@ date version nhelpfiles nobjects ndatasets Rlines srclines
 "2024-11-19" "3.3-3" 269 726 0 39066 1155
 "2025-01-20" "3.3-4" 269 726 0 39070 1155
 "2025-03-22" "3.3-5" 269 726 0 39076 1155
-"2025-04-18" "3.3-5.002" 269 727 0 39131 1155
+"2025-04-19" "3.3-5.002" 269 729 0 39061 1155

inst/info/packagesizes.txt

@@ -19,4 +19,4 @@ date version nhelpfiles nobjects ndatasets Rlines srclines
 "2024-11-19" "3.3-3" 269 726 0 39066 1155
 "2025-01-20" "3.3-4" 269 726 0 39070 1155
 "2025-03-22" "3.3-5" 269 726 0 39076 1155
-"2025-04-18" "3.3-5.002" 269 727 0 39131 1155
+"2025-04-19" "3.3-5.002" 269 729 0 39061 1155

man/spatstat.model-internal.Rd

@@ -17,7 +17,6 @@
 \alias{coef.summary.ppm}
 \alias{coef.summary.slrm}
 \alias{coef.vblogit} 
-\alias{condSimCox}
 \alias{damaged.ppm}
 \alias{data.mppm}
 \alias{deltasuffstat}
@@ -94,7 +93,9 @@
 \alias{is.multitype.slrm}    
 \alias{is.poisson.mppm}
 \alias{is.poisson.rppm}
+\alias{is.poissonclusterprocess.clusterprocess}
 \alias{Kpcf.kppm}
+\alias{Kmodel.clusterprocess}
 \alias{Kmodel.slrm}
 \alias{killinteraction}
 \alias{kppmComLik}
@@ -129,6 +130,7 @@
 \alias{parameters.detpointprocfamily}
 \alias{PairPotentialType}
 \alias{partialModelMatrix}
+\alias{pcfmodel.clusterprocess}
 \alias{pcfmodel.slrm}
 \alias{ploterodewin}
 \alias{ploterodeimage}
@@ -260,8 +262,6 @@ check.separable(dmat, covname, isconstant, fatal)
 \method{coef}{summary.ppm}(object, \dots)
 \method{coef}{summary.slrm}(object, \dots)
 \method{coef}{vblogit}(object, \dots) 
-condSimCox(object, nsim, \dots, window, n.cond, w.cond,
-           giveup, maxchunk, saveLambda, verbose, drop)
 damaged.ppm(object)
 data.mppm(x)
 deltasuffstat(model, \dots,
@@ -346,7 +346,9 @@ is.mppm(x)
 \method{is.multitype}{slrm}(X, \dots)
 \method{is.poisson}{mppm}(x)
 \method{is.poisson}{rppm}(x)
+\method{is.poissonclusterprocess}{clusterprocess}(model)
 Kpcf.kppm(model, what)
+\method{Kmodel}{clusterprocess}(model, \dots)
 \method{Kmodel}{slrm}(model, \dots)
 killinteraction(model)
 kppmComLik(X, Xname, po, clusters, control, stabilize, weightfun, rmax,
@@ -419,6 +421,7 @@ oversize.quad(Q, \dots, nU, nX, p)
 PairPotentialType(pairpot)
 \method{parameters}{detpointprocfamily}(model, \dots)
 partialModelMatrix(X,D,model,callstring,\dots)
+\method{pcfmodel}{clusterprocess}(model, \dots)
 \method{pcfmodel}{slrm}(model, \dots)
 ploterodewin(W1, W2, col.edge, col.inside, do.plot, \dots)
 ploterodeimage(W, Z, \dots, Wcol, rangeZ, colsZ, do.plot)