rmh.ppm: Simulate from a Fitted Point Process Model

Description

Given a point process model fitted to data, generate a random simulation of the model, using the Metropolis-Hastings algorithm.

Usage

## S3 method for class 'ppm':
rmh(model, start=NULL,
                    control=default.rmhcontrol(model),
                    ...,
                    project=TRUE, verbose=TRUE)

Arguments

model

A fitted point process model (object of class "ppm", see ppm.object) which it is desired to simulate. This fitted model is usually the result of a call to

start

Data determining the initial state of the Metropolis-Hastings algorithm. See rmhstart for description of these arguments. Defaults to list(x.start=data.ppm(model))

control

Data controlling the iterative behaviour of the Metropolis-Hastings algorithm. See rmhcontrol for description of these arguments.

...

Further arguments passed to rmhcontrol, or to rmh.default, or to covariate functions in the model.

project

Logical flag indicating what to do if the fitted model is invalid (in the sense that the values of the fitted coefficients do not specify a valid point process). If project=TRUE the closest valid model will be simulated; if

verbose

Logical flag indicating whether to print progress reports.

Value

A point pattern (an object of class "ppp"; see ppp.object).

Warnings

See Warnings in rmh.default.

Details

This function generates simulated realisations from a point process model that has been fitted to point pattern data. It is a method for the generic function rmh for the class "ppm" of fitted point process models. To simulate other kinds of point process models, see rmh or rmh.default.

The argument model describes the fitted model. It must be an object of class "ppm" (see ppm.object), and will typically be the result of a call to the point process model fitting function ppm.

The current implementation enables simulation from any fitted model involving the interactions DiggleGratton, Geyer, MultiStrauss, MultiStraussHard, PairPiece, Poisson, Strauss, StraussHard and Softcore, including nonstationary models. See the examples.

It is possible that the fitted coefficients of a point process model may be ``illegal'', i.e. that there may not exist a mathematically well-defined point process with the given parameter values. For example, a Strauss process with interaction parameter $\gamma > 1$ does not exist, but the model-fitting procedure used in ppm will sometimes produce values of $\gamma$ greater than 1. In such cases, if project=FALSE then an error will occur, while if project=TRUE then rmh.ppm will find the nearest legal model and simulate this model instead. (The nearest legal model is obtained by projecting the vector of coefficients onto the set of valid coefficient vectors. The result is usually the Poisson process with the same fitted intensity.) The arguments start and control are lists of parameters determining the initial state and the iterative behaviour, respectively, of the Metropolis-Hastings algorithm.

The argument start is passed directly to rmhstart. See rmhstart for details of the parameters of the initial state, and their default values.

The argument control is first passed to rmhcontrol. Then if any additional arguments ... are given, update.rmhcontrol is called to update the parameter values. See rmhcontrol for details of the iterative behaviour parameters, and default.rmhcontrol for their default values.

Note that if you specify expansion of the simulation window using the parameter expand (so that the model will be simulated on a window larger than the original data window) then the model must be capable of extrapolation to this larger window. This is usually not possible for models which depend on external covariates, because the domain of a covariate image is usually the same as the domain of the fitted model. After extracting the relevant information from the fitted model object model, rmh.ppm invokes the default rmh algorithm rmh.default, unless the model is Poisson. If the model is Poisson then the Metropolis-Hastings algorithm is not needed, and the model is simulated directly, using one of rpoispp, rmpoispp, rpoint or rmpoint.

See rmh.default for further information about the implementation, or about the Metropolis-Hastings algorithm.

Examples

Run this code

live <- interactive()
   op <- spatstat.options()
   spatstat.options(rmh.nrep=1e5)
   Nrep <- 1e5

   <testonly>spatstat.options(rmh.nrep=10, npixel=10)
     Nrep <- 10</testonly>

   data(swedishpines)
   X <- swedishpines
   plot(X, main="Swedish Pines data")

   # Poisson process
   fit <- ppm(X, ~1, Poisson())
   Xsim <- rmh(fit)
   if(live) plot(Xsim, main="simulation from fitted Poisson model")

   # Strauss process   
   fit <- ppm(X, ~1, Strauss(r=7))
   Xsim <- rmh(fit)
   if(live) plot(Xsim, main="simulation from fitted Strauss model")

   # Strauss process simulated on a larger window
   # then clipped to original window
   Xsim <- rmh(fit, control=list(nrep=Nrep, expand=2, periodic=TRUE))
   Xsim <- rmh(fit, nrep=Nrep, expand=2, periodic=TRUE)
   
   # Strauss - hard core process
   fit <- ppm(X, ~1, StraussHard(r=7,hc=2))
   Xsim <- rmh(fit, start=list(n.start=X$n))
   if(live) plot(Xsim, main="simulation from fitted Strauss hard core model")

   # Geyer saturation process
   fit <- ppm(X, ~1, Geyer(r=7,sat=2))
   Xsim <- rmh(fit, start=list(n.start=X$n))
   if(live) plot(Xsim, main="simulation from fitted Geyer model")

   # Area-interaction process
   fit <- ppm(X, ~1, AreaInter(r=7))
   Xsim <- rmh(fit, start=list(n.start=X$n))
   if(live) plot(Xsim, main="simulation from fitted area-interaction model")

   # soft core interaction process
   Q <- quadscheme(X, nd=50)
   fit <- ppm(Q, ~1, Softcore(kappa=0.1), correction="isotropic")
   Xsim <- rmh(fit, start=list(n.start=X$n))
   if(live) plot(Xsim, main="simulation from fitted Soft Core model")

   data(cells)
   if(live) plot(cells)
   # Diggle-Gratton pairwise interaction model
   fit <- ppm(cells, ~1, DiggleGratton(0.05, 0.1))
   Xsim <- rmh(fit, start=list(n.start=cells$n))
   if(live) plot(Xsim, main="simulation from fitted Diggle-Gratton model")

   X <- rSSI(0.05, 100)
   if(live) plot(X, main="new data")

   # piecewise-constant pairwise interaction function
   fit <- ppm(X, ~1, PairPiece(seq(0.02, 0.1, by=0.01)))
   Xsim <- rmh(fit)
   if(live) plot(Xsim, main="simulation from fitted pairwise model")

   # marked point pattern
   data(amacrine)
   Y <- amacrine
   if(live) plot(Y, main="Amacrine data")

   # marked Poisson models 
   fit <- ppm(Y)
   Ysim <- rmh(fit)
   if(live) plot(Ysim, main="simulation from ppm(Y)")

   fit <- ppm(Y,~marks)
   Ysim <- rmh(fit)
   if(live) plot(Ysim, main="simulation from ppm(Y, ~marks)")

   fit <- ppm(Y,~polynom(x,y,2))
   Ysim <- rmh(fit)
   if(live) plot(Ysim, main="simulation from ppm(Y, ~polynom(x,y,2))")

   fit <- ppm(Y,~marks+polynom(x,y,2))
   Ysim <- rmh(fit)
   if(live) plot(Ysim, main="simulation from ppm(Y, ~marks+polynom(x,y,2))")

   fit <- ppm(Y,~marks*polynom(x,y,2))
   Ysim <- rmh(fit)
   if(live) plot(Ysim, main="simulation from ppm(Y, ~marks*polynom(x,y,2))")

   # multitype Strauss models
   MS <- MultiStrauss(types = levels(Y$marks),
                      radii=matrix(0.07, ncol=2, nrow=2))
   fit <- ppm(Y, ~marks, MS)
   Ysim <- rmh(fit)
   if(live) plot(Ysim, main="simulation from fitted Multitype Strauss")

   fit <- ppm(Y,~marks*polynom(x,y,2), MS)
   Ysim <- rmh(fit)
   if(live) plot(Ysim, main="simulation from fitted inhomogeneous Multitype Strauss")

   spatstat.options(op)

Run the code above in your browser using DataLab