Package 'spatstat.random'

The spatstat.random Package

Description

The spatstat.random package belongs to the spatstat family of packages. It contains the functionality for generating random spatial patterns and simulation of random point processes.

Details

spatstat is a family of R packages for the statistical analysis of spatial data. Its main focus is the analysis of spatial patterns of points in two-dimensional space.

This sub-package spatstat.random contains the functions that perform random generation of spatial patterns and simulation of random point processes:

• generation of random spatial patterns of points according to many simple rules (completely random patterns, random grids, systematic random points);

• randomised alteration of spatial patterns (thinning, random shifting, jittering, random labelling);

• generation of quasirandom patterns;

• direct simulation of random point processes (Poisson process, binomial process, cell process, simple sequential inhibition, Matern inhibition models, log-Gaussian Cox processes;

• simulation of Neyman-Scott cluster processes (truncated direct algorithm, Brix-Kendall and hybrid algorithms) and product shot noise cluster processes;

• simulation of Gibbs point processes (Metropolis-Hastings birth-death-shift algorithm, perfect simulation/ dominated coupling from the past, alternating Gibbs sampler).

Some other types of spatial object are also supported:

• generation of random patterns of points in 3 dimensions;

• generation of random spatial patterns of line segments;

• generation of random tessellations;

• generation of random images (random noise, random mosaics).

(Functions for linear networks are in the separate sub-package spatstat.linnet.)

Structure of the spatstat family

The spatstat family of packages is designed to support a complete statistical analysis of spatial data. It supports

• creation, manipulation and plotting of point patterns;

• exploratory data analysis;

• spatial random sampling;

• simulation of point process models;

• parametric model-fitting;

• non-parametric smoothing and regression;

• formal inference (hypothesis tests, confidence intervals);

• model diagnostics.

The orginal spatstat package grew to be very large. It has now been divided into several sub-packages:

• spatstat.utils containing basic utilities

• spatstat.sparse containing linear algebra utilities

• spatstat.data containing datasets

• spatstat.univar containing functions for estimating probability distributions of random variables

• spatstat.geom containing geometrical objects and geometrical operations

• spatstat.random containing functionality for simulation and random generation

• spatstat.explore containing the main functionality for exploratory data analysis and nonparametric analysis

• spatstat.model containing the main functionality for parametric modelling and formal inference for spatial data

• spatstat.linnet containing functions for spatial data on a linear network

• spatstat, which simply loads the other sub-packages listed above, and provides documentation.

When you install spatstat, these sub-packages are also installed. Then if you load the spatstat package by typing library(spatstat), the other sub-packages listed above will automatically be loaded or imported.

For an overview of all the functions available in the sub-packages of spatstat, see the help file for "spatstat-package" in the spatstat package.

Additionally there are several extension packages:

• spatstat.gui for interactive graphics

• spatstat.local for local likelihood (including geographically weighted regression)

• spatstat.Knet for additional, computationally efficient code for linear networks

• spatstat.sphere (under development) for spatial data on a sphere, including spatial data on the earth's surface

The extension packages must be installed separately and loaded explicitly if needed. They also have separate documentation.

Functionality in spatstat.random

Following is a list of the functionality provided in the spatstat.random package only.

To simulate a random point pattern:

runifpoint	generate $n$ independent uniform random points
rpoint	generate $n$ independent random points
rmpoint	generate $n$ independent multitype random points
rpoispp	simulate the (in)homogeneous Poisson point process
rmpoispp	simulate the (in)homogeneous multitype Poisson point process
runifdisc	generate $n$ independent uniform random points in disc
rstrat	stratified random sample of points
rMaternI	simulate the Matern Model I inhibition process
rMaternII	simulate the Matern Model II inhibition process
rSSI	simulate Simple Sequential Inhibition process
rStrauss	simulate Strauss process (perfect simulation)
rHardcore	simulate Hard Core process (perfect simulation)
rStraussHard	simulate Strauss-hard core process (perfect simulation)
rDiggleGratton	simulate Diggle-Gratton process (perfect simulation)
rDGS	simulate Diggle-Gates-Stibbard process (perfect simulation)
rPenttinen	simulate Penttinen process (perfect simulation)
rNeymanScott	simulate a general Neyman-Scott process
rPoissonCluster	simulate a general Poisson cluster process
rMatClust	simulate the Matern Cluster process
rThomas	simulate the Thomas process
rGaussPoisson	simulate the Gauss-Poisson cluster process
rCauchy	simulate Neyman-Scott Cauchy cluster process
rVarGamma	simulate Neyman-Scott Variance Gamma cluster process
rthin	random thinning
rcell	simulate the Baddeley-Silverman cell process
rmh	simulate Gibbs point process using Metropolis-Hastings
runifpointOnLines	generate $n$ random points along specified line segments
rpoisppOnLines	generate Poisson random points along specified line segments

To randomly change an existing point pattern:

rshift	random shifting of points
rthin	random thinning
rlabel	random (re)labelling of a multitype point pattern
quadratresample	block resampling

See also rjitter and rexplode in the spatstat.geom package.

Random pixel images:

An object of class "im" represents a pixel image.

rnoise

random pixel noise

Line segment patterns

An object of class "psp" represents a pattern of straight line segments.

rpoisline

generate a realisation of the Poisson line process inside a window

Tessellations

An object of class "tess" represents a tessellation.

rpoislinetess

generate tessellation using Poisson line process

Three-dimensional point patterns

An object of class "pp3" represents a three-dimensional point pattern in a rectangular box. The box is represented by an object of class "box3".

runifpoint3	generate uniform random points in 3-D
rpoispp3	generate Poisson random points in 3-D

Multi-dimensional space-time point patterns

An object of class "ppx" represents a point pattern in multi-dimensional space and/or time.

runifpointx	generate uniform random points
rpoisppx	generate Poisson random points

Probability Distributions

rknn	theoretical distribution of nearest neighbour distance
dmixpois	mixed Poisson distribution

Simulation

There are many ways to generate a random point pattern, line segment pattern, pixel image or tessellation in spatstat.

Random point patterns:

runifpoint	generate $n$ independent uniform random points
rpoint	generate $n$ independent random points
rmpoint	generate $n$ independent multitype random points
rpoispp	simulate the (in)homogeneous Poisson point process
rmpoispp	simulate the (in)homogeneous multitype Poisson point process
runifdisc	generate $n$ independent uniform random points in disc
rstrat	stratified random sample of points
rMaternI	simulate the Matern Model I inhibition process
rMaternII	simulate the Matern Model II inhibition process
rSSI	simulate Simple Sequential Inhibition process
rHardcore	simulate hard core process (perfect simulation)
rStrauss	simulate Strauss process (perfect simulation)
rStraussHard	simulate Strauss-hard core process (perfect simulation)
rDiggleGratton	simulate Diggle-Gratton process (perfect simulation)
rDGS	simulate Diggle-Gates-Stibbard process (perfect simulation)
rPenttinen	simulate Penttinen process (perfect simulation)
rNeymanScott	simulate a general Neyman-Scott process
rMatClust	simulate the Matern Cluster process
rThomas	simulate the Thomas process
rLGCP	simulate the log-Gaussian Cox process
rGaussPoisson	simulate the Gauss-Poisson cluster process
rCauchy	simulate Neyman-Scott process with Cauchy clusters
rVarGamma	simulate Neyman-Scott process with Variance Gamma clusters
rcell	simulate the Baddeley-Silverman cell process
runifpointOnLines	generate $n$ random points along specified line segments
rpoisppOnLines	generate Poisson random points along specified line segments

Resampling a point pattern:

quadratresample	block resampling
rshift	random shifting of (subsets of) points
rthin	random thinning

Other random patterns:

rpoisline	simulate the Poisson line process within a window
rpoislinetess	generate random tessellation using Poisson line process
rMosaicSet	generate random set by selecting some tiles of a tessellation
rMosaicField	generate random pixel image by assigning random values in each tile of a tessellation

Resampling and randomisation procedures

You can build your own tests based on randomisation and resampling using the following capabilities:

quadratresample	block resampling
rshift	random shifting of (subsets of) points
rthin	random thinning

Licence

This library and its documentation are usable under the terms of the "GNU General Public License", a copy of which is distributed with the package.

Acknowledgements

Kasper Klitgaard Berthelsen, Ya-Mei Chang, Tilman Davies, Ute Hahn, Abdollah Jalilian, Dominic Schuhmacher and Rasmus Waagepetersen made substantial contributions of code.

For comments, corrections, bug alerts and suggestions, we thank Monsuru Adepeju, Corey Anderson, Ang Qi Wei, Ryan Arellano, Jens Astrom, Robert Aue, Marcel Austenfeld, Sandro Azaele, Malissa Baddeley, Guy Bayegnak, Colin Beale, Melanie Bell, Thomas Bendtsen, Ricardo Bernhardt, Andrew Bevan, Brad Biggerstaff, Anders Bilgrau, Leanne Bischof, Christophe Biscio, Roger Bivand, Jose M. Blanco Moreno, Florent Bonneu, Jordan Brown, Ian Buller, Julian Burgos, Simon Byers, Ya-Mei Chang, Jianbao Chen, Igor Chernayavsky, Y.C. Chin, Bjarke Christensen, Lucia Cobo Sanchez, Jean-Francois Coeurjolly, Kim Colyvas, Hadrien Commenges, Rochelle Constantine, Robin Corria Ainslie, Richard Cotton, Marcelino de la Cruz, Peter Dalgaard, Mario D'Antuono, Sourav Das, Peter Diggle, Patrick Donnelly, Ian Dryden, Stephen Eglen, Ahmed El-Gabbas, Belarmain Fandohan, Olivier Flores, David Ford, Peter Forbes, Shane Frank, Janet Franklin, Funwi-Gabga Neba, Oscar Garcia, Agnes Gault, Jonas Geldmann, Marc Genton, Shaaban Ghalandarayeshi, Jason Goldstick, Pavel Grabarnik, C. Graf, Ute Hahn, Andrew Hardegen, Martin Bogsted Hansen, Martin Hazelton, Juha Heikkinen, Mandy Hering, Markus Herrmann, Maximilian Hesselbarth, Paul Hewson, Hamidreza Heydarian, Kurt Hornik, Philipp Hunziker, Jack Hywood, Ross Ihaka, Cenk Icos, Aruna Jammalamadaka, Robert John-Chandran, Devin Johnson, Mahdieh Khanmohammadi, Bob Klaver, Lily Kozmian-Ledward, Peter Kovesi, Mike Kuhn, Jeff Laake, Robert Lamb, Frederic Lavancier, Tom Lawrence, Tomas Lazauskas, Jonathan Lee, George Leser, Angela Li, Li Haitao, George Limitsios, Andrew Lister, Nestor Luambua, Ben Madin, Martin Maechler, Kiran Marchikanti, Jeff Marcus, Robert Mark, Peter McCullagh, Monia Mahling, Jorge Mateu Mahiques, Ulf Mehlig, Frederico Mestre, Sebastian Wastl Meyer, Mi Xiangcheng, Lore De Middeleer, Robin Milne, Enrique Miranda, Jesper Moller, Annie Mollie, Ines Moncada, Mehdi Moradi, Virginia Morera Pujol, Erika Mudrak, Gopalan Nair, Nader Najari, Nicoletta Nava, Linda Stougaard Nielsen, Felipe Nunes, Jens Randel Nyengaard, Jens Oehlschlaegel, Thierry Onkelinx, Sean O'Riordan, Evgeni Parilov, Jeff Picka, Nicolas Picard, Tim Pollington, Mike Porter, Sergiy Protsiv, Adrian Raftery, Ben Ramage, Pablo Ramon, Xavier Raynaud, Nicholas Read, Matt Reiter, Ian Renner, Tom Richardson, Brian Ripley, Ted Rosenbaum, Barry Rowlingson, Jason Rudokas, Tyler Rudolph, John Rudge, Christopher Ryan, Farzaneh Safavimanesh, Aila Sarkka, Cody Schank, Katja Schladitz, Sebastian Schutte, Bryan Scott, Olivia Semboli, Francois Semecurbe, Vadim Shcherbakov, Shen Guochun, Shi Peijian, Harold-Jeffrey Ship, Tammy L Silva, Ida-Maria Sintorn, Yong Song, Malte Spiess, Mark Stevenson, Kaspar Stucki, Jan Sulavik, Michael Sumner, P. Surovy, Ben Taylor, Thordis Linda Thorarinsdottir, Leigh Torres, Berwin Turlach, Torben Tvedebrink, Kevin Ummer, Medha Uppala, Andrew van Burgel, Tobias Verbeke, Mikko Vihtakari, Alexendre Villers, Fabrice Vinatier, Maximilian Vogtland, Sasha Voss, Sven Wagner, Hao Wang, H. Wendrock, Jan Wild, Carl G. Witthoft, Selene Wong, Maxime Woringer, Luke Yates, Mike Zamboni and Achim Zeileis.

Author(s)

Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

---

Convert Data To Class owin

Description

Converts data specifying an observation window in any of several formats, into an object of class "owin".

Usage

## S3 method for class 'rmhmodel'
as.owin(W, ..., fatal=FALSE)

Arguments

W	Data specifying an observation window, in any of several formats described under Details below.
fatal	Logical value determining what to do if the data cannot be converted to an observation window. See Details.
...	Ignored.

Details

The class "owin" is a way of specifying the observation window for a point pattern. See owin.object for an overview.

The generic function as.owin converts data in any of several formats into an object of class "owin" for use by the spatstat package. The function as.owin is generic, with methods for different classes of objects, and a default method.

The argument W may be

• an object of class "owin"

• a structure with entries xrange, yrange specifying the $x$ and $y$ dimensions of a rectangle

• a structure with entries named xmin, xmax, ymin, ymax (in any order) specifying the $x$ and $y$ dimensions of a rectangle. This will accept objects of class bbox in the sf package.

• a numeric vector of length 4 (interpreted as (xmin, xmax, ymin, ymax) in that order) specifying the $x$ and $y$ dimensions of a rectangle

• a structure with entries named xl, xu, yl, yu (in any order) specifying the $x$ and $y$ dimensions of a rectangle as (xmin, xmax) = (xl, xu) and (ymin, ymax) = (yl, yu). This will accept objects of class spp used in the Venables and Ripley spatial package.

• an object of class "ppp" representing a point pattern. In this case, the object's window structure will be extracted.

• an object of class "psp" representing a line segment pattern. In this case, the object's window structure will be extracted.

• an object of class "tess" representing a tessellation. In this case, the object's window structure will be extracted.

• an object of class "quad" representing a quadrature scheme. In this case, the window of the data component will be extracted.

• an object of class "im" representing a pixel image. In this case, a window of type "mask" will be returned, with the same pixel raster coordinates as the image. An image pixel value of NA, signifying that the pixel lies outside the window, is transformed into the logical value FALSE, which is the corresponding convention for window masks.

• an object of class "ppm", "kppm", "slrm" or "dppm" representing a fitted point process model. In this case, if from="data" (the default), as.owin extracts the original point pattern data to which the model was fitted, and returns the observation window of this point pattern. If from="covariates" then as.owin extracts the covariate images to which the model was fitted, and returns a binary mask window that specifies the pixel locations.

• an object of class "lpp" representing a point pattern on a linear network. In this case, as.owin extracts the linear network and returns a window containing this network.

• an object of class "lppm" representing a fitted point process model on a linear network. In this case, as.owin extracts the linear network and returns a window containing this network.

• A data.frame with exactly three columns. Each row of the data frame corresponds to one pixel. Each row contains the $x$ and $y$ coordinates of a pixel, and a logical value indicating whether the pixel lies inside the window.

• A data.frame with exactly two columns. Each row of the data frame contains the $x$ and $y$ coordinates of a pixel that lies inside the window.

• an object of class "distfun", "nnfun" or "funxy" representing a function of spatial location, defined on a spatial domain. The spatial domain of the function will be extracted.

• an object of class "rmhmodel" representing a point process model that can be simulated using rmh. The window (spatial domain) of the model will be extracted. The window may be NULL in some circumstances (indicating that the simulation window has not yet been determined). This is not treated as an error, because the argument fatal defaults to FALSE for this method.

• an object of class "layered" representing a list of spatial objects. See layered. In this case, as.owin will be applied to each of the objects in the list, and the union of these windows will be returned.

• an object of another suitable class from another package. For full details, see vignette('shapefiles').

If the argument W is not in one of these formats and cannot be converted to a window, then an error will be generated (if fatal=TRUE) or a value of NULL will be returned (if fatal=FALSE).

When W is a data frame, the argument step can be used to specify the pixel grid spacing; otherwise, the spacing will be guessed from the data.

Value

An object of class "owin" (see owin.object) specifying an observation window.

Author(s)

Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

See Also

as.owin, as.owin.ppm, as.owin.lpp.

owin.object, owin.

Additional methods for as.owin may be provided by other packages outside the spatstat family.

Examples

m <- rmhmodel(cif='poisson', par=list(beta=1), w=square(2))
  as.owin(m)

---

Field of clusters

Description

Calculate the superposition of cluster kernels at the location of a point pattern.

Usage

clusterfield(model, locations = NULL, ...)

  ## S3 method for class 'character'
clusterfield(model, locations = NULL, ...)

  ## S3 method for class 'function'
clusterfield(model, locations = NULL, ..., mu = NULL)

Arguments

model	Cluster model. Either a fitted cluster model (object of class "kppm"), a character string specifying the type of cluster model, or a function defining the cluster kernel. See Details.
locations	A point pattern giving the locations of the kernels. Defaults to the centroid of the observation window for the "kppm" method and to the center of a unit square otherwise.
...	Additional arguments passed to density.ppp or the cluster kernel. See Details.
mu	Mean number of offspring per cluster. A single number or a pixel image.

Details

The function clusterfield is generic, with methods for "character" and "function" (described here) and a method for "kppm" (described in clusterfield.kppm).

The calculations are performed by density.ppp and ... arguments are passed thereto for control over the pixel resolution etc. (These arguments are then passed on to pixellate.ppp and as.mask.)

For the method clusterfield.function, the given kernel function should accept vectors of x and y coordinates as its first two arguments. Any additional arguments may be passed through the ....

The method clusterfield.function also accepts the optional parameter mu (defaulting to 1) specifying the mean number of points per cluster (as a numeric) or the inhomogeneous reference cluster intensity (as an "im" object or a function(x,y)). The interpretation of mu is as explained in the simulation functions referenced in the See Also section below.

For the method clusterfield.character, the argument model must be one of the following character strings: model="Thomas" for the Thomas process, model="MatClust" for the Matern cluster process, model="Cauchy" for the Neyman-Scott cluster process with Cauchy kernel, or model="VarGamma" for the Neyman-Scott cluster process with Variance Gamma kernel. For all these models the parameter scale is required and passed through ... as well as the parameter nu when model="VarGamma". This method calls clusterfield.function so the parameter mu may also be passed through ... and will be interpreted as explained above.

Value

A pixel image (object of class "im").

Author(s)

Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

See Also

clusterfield.kppm.

density.ppp and kppm.

Simulation algorithms for cluster models: rCauchy rMatClust rThomas rVarGamma

Examples

# method for functions
  kernel <- function(x,y,scal) { 
      r <- sqrt(x^2 + y^2)
      ifelse(r > 0,
             dgamma(r, shape=5, scale=scal)/(2 * pi * r),
             0)               
  }
  X <- runifpoint(10)
  clusterfield(kernel, X, scal=0.05)

---

Extract Cluster Offspring Kernel

Description

Given a cluster point process model, this command returns the probability density of the cluster offspring.

Usage

clusterkernel(model, ...)
## S3 method for class 'character'
clusterkernel(model, ...)

Arguments

model	Cluster model. Either a fitted cluster or Cox model (object of class "kppm"), or a character string specifying the type of cluster model.
...	Parameter values for the model, when model is a character string.

Details

Given a specification of a cluster point process model, this command returns a function(x,y) giving the two-dimensional probability density of the cluster offspring points assuming a cluster parent located at the origin.

The function clusterkernel is generic, with methods for class "character" (described here) and "kppm" (described in clusterkernel.kppm).

Value

A function in the R language with arguments x,y,....

Author(s)

Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

See Also

clusterkernel.kppm), clusterfield, kppm.

Examples

f <- clusterkernel("Thomas", kappa=10, scale=0.5)
  f(0.1, 0.2)

---

Compute or Extract Effective Range of Cluster Kernel

Description

Given a cluster point process model, this command returns a value beyond which the the probability density of the cluster offspring is neglible.

Usage

clusterradius(model, ...)

## S3 method for class 'character'
clusterradius(model, ..., thresh = NULL, precision = FALSE)

Arguments

model	Cluster model. Either a fitted cluster or Cox model (object of class "kppm"), or a character string specifying the type of cluster model.
...	Parameter values for the model, when model is a character string.
thresh	Numerical threshold relative to the cluster kernel value at the origin (parent location) determining when the cluster kernel will be considered neglible. A sensible default is provided.
precision	Logical. If precision=TRUE the precision of the calculated range is returned as an attribute to the range. See details.

Details

Given a cluster model this function by default returns the effective range of the model with the given parameters as used in spatstat. For the Matern cluster model (see e.g. rMatClust) this is simply the finite radius of the offsring density given by the paramter scale irrespective of other options given to this function. The remaining models in spatstat have infinite theoretical range, and an effective finite value is given as follows: For the Thomas model (see e.g. rThomas the default is 4*scale where scale is the scale or standard deviation parameter of the model. If thresh is given the value is instead found as described for the other models below.

For the Cauchy model (see e.g. rCauchy) and the Variance Gamma (Bessel) model (see e.g. rVarGamma) the value of thresh defaults to 0.001, and then this is used to compute the range numerically as follows. If $k(x,y)=k_0(r)$ with $r=\sqrt(x^2+y^2)$ denotes the isotropic cluster kernel then $f(r) = 2 \pi r k_0(r)$ is the density function of the offspring distance from the parent. The range is determined as the value of $r$ where $f(r)$ falls below thresh times $k_0(r)$.

If precision=TRUE the precision related to the chosen range is returned as an attribute. Here the precision is defined as the polar integral of the kernel from distance 0 to the calculated range. Ideally this should be close to the value 1 which would be obtained for the true theretical infinite range.

Value

A positive numeric.

Additionally, the precision related to this range value is returned as an attribute "prec", if precision=TRUE.

Author(s)

Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

See Also

clusterkernel, kppm, rMatClust, rThomas, rCauchy, rVarGamma, rNeymanScott.

Examples

clusterradius("Thomas", scale = .1)
  clusterradius("Thomas", scale = .1, thresh = 0.001)
  clusterradius("VarGamma", scale = .1, nu = 2, precision = TRUE)

---

Default Expansion Rule for Simulation of Model

Description

Defines the default expansion window or expansion rule for simulation of a point process model.

Usage

default.expand(object, m=2, epsilon=1e-6, w=Window(object))

Arguments

object	A point process model (object of class "ppm" or "rmhmodel").
m	A single numeric value. The window will be expanded by a distance m * reach(object) along each side.
epsilon	Threshold argument passed to reach to determine reach(object).
w	Optional. The un-expanded window in which the model is defined. The resulting simulated point patterns will lie in this window.

Details

This function computes a default value for the expansion rule (the argument expand in rmhcontrol) given a fitted point process model object. This default is used by rmh, simulate.ppm, envelope, qqplot.ppm, and other functions.

Suppose we wish to generate simulated realisations of a fitted point process model inside a window w. It is advisable to first simulate the pattern on a larger window, and then clip it to the original window w. This avoids edge effects in the simulation. It is called expansion of the simulation window.

Accordingly, for the Metropolis-Hastings simulation algorithm rmh, the algorithm control parameters specified by rmhcontrol include an argument expand that determines the expansion of the simulation window.

The function default.expand determines the default expansion rule for a fitted point process model object.

If the model is Poisson, then no expansion is necessary. No expansion is performed by default, and default.expand returns a rule representing no expansion. The simulation window is the original window w = Window(object).

If the model depends on external covariates (i.e.\ covariates other than the Cartesian covariates x and y and the marks) then no expansion is feasible, in general, because the spatial domain of the covariates is not guaranteed to be large enough. default.expand returns a rule representing no expansion. The simulation window is the original window w = Window(object).

If the model depends on the Cartesian covariates x and y, it would be feasible to expand the simulation window, and this was the default for spatstat version 1.24-1 and earlier. However this sometimes produces artefacts (such as an empty point pattern) or memory overflow, because the fitted trend, extrapolated outside the original window of the data, may become very large. In spatstat version 1.24-2 and later, the default rule is not to expand if the model depends on x or y. Again default.expand returns a rule representing no expansion.

Otherwise, expansion will occur. The original window w = Window(object) is expanded by a distance m * rr, where rr is the interaction range of the model, computed by reach. If w is a rectangle then each edge of w is displaced outward by distance m * rr. If w is not a rectangle then w is dilated by distance m * rr using dilation.

Value

A window expansion rule (object of class "rmhexpand").

Author(s)

Adrian Baddeley [email protected] and Rolf Turner [email protected]

See Also

rmhexpand, rmhcontrol, rmh, envelope, qqplot.ppm

Examples

if(require(spatstat.model)) {
  fit <- ppm(cells ~1, Strauss(0.07))
  default.expand(fit)
}
  mod <- rmhmodel(cif="strauss", par=list(beta=100, gamma=0.5, r=0.07))
  default.expand(mod)

---

Set Default Control Parameters for Metropolis-Hastings Algorithm.

Description

For a Gibbs point process model (either a fitted model, or a model specified by its parameters), this command sets appropriate default values of the parameters controlling the iterative behaviour of the Metropolis-Hastings algorithm.

Usage

default.rmhcontrol(model, w=NULL)

Arguments

model	A fitted point process model (object of class "ppm") or a description of a Gibbs point process model (object of class "rmhmodel").
w	Optional. Window for the resulting simulated patterns.

Details

This function sets the values of the parameters controlling the iterative behaviour of the Metropolis-Hastings simulation algorithm. It uses default values that would be appropriate for the fitted point process model model.

The expansion parameter expand is set to default.expand(model, w).

All other parameters revert to their defaults given in rmhcontrol.default.

See rmhcontrol for the full list of control parameters. To override default parameters, use update.rmhcontrol.

Value

An object of class "rmhcontrol". See rmhcontrol.

Author(s)

Adrian Baddeley [email protected] and Rolf Turner [email protected]

See Also

rmhcontrol, update.rmhcontrol, ppm, default.expand

Examples

if(require(spatstat.model)) {
  fit <- ppm(cells, ~1, Strauss(0.1))
  default.rmhcontrol(fit)
  default.rmhcontrol(fit, w=square(2))
}
   m <- rmhmodel(cif='strauss',
                 par=list(beta=100, gamma=0.5, r=0.1),
                 w=unit.square())
  default.rmhcontrol(m)
  default.rmhcontrol(m, w=square(2))

---

Mixed Poisson Distribution

Description

Density, distribution function, quantile function and random generation for a mixture of Poisson distributions.

Usage

dmixpois(x, mu, sd, invlink = exp, GHorder = 5)
pmixpois(q, mu, sd, invlink = exp, lower.tail = TRUE, GHorder = 5)
qmixpois(p, mu, sd, invlink = exp, lower.tail = TRUE, GHorder = 5)
rmixpois(n, mu, sd, invlink = exp)

Arguments

x	vector of (non-negative integer) quantiles.
q	vector of quantiles.
p	vector of probabilities.
n	number of random values to return.
mu	Mean of the linear predictor. A single numeric value.
sd	Standard deviation of the linear predictor. A single numeric value.
invlink	Inverse link function. A function in the R language, used to transform the linear predictor into the parameter lambda of the Poisson distribution.
lower.tail	Logical. If TRUE (the default), probabilities are $P[X \le x]$, otherwise, $P[X > x]$.
GHorder	Number of quadrature points in the Gauss-Hermite quadrature approximation. A small positive integer.

Details

These functions are analogous to dpois ppois, qpois and rpois except that they apply to a mixture of Poisson distributions.

In effect, the Poisson mean parameter lambda is randomised by setting lambda = invlink(Z) where Z has a Gaussian $N(\mu,\sigma^2)$ distribution. The default is invlink=exp which means that lambda is lognormal. Set invlink=I to assume that lambda is approximately Normal.

For dmixpois, pmixpois and qmixpois, the probability distribution is approximated using Gauss-Hermite quadrature. For rmixpois, the deviates are simulated exactly.

Value

Numeric vector: dmixpois gives probability masses, ppois gives cumulative probabilities, qpois gives (non-negative integer) quantiles, and rpois generates (non-negative integer) random deviates.

Author(s)

Adrian Baddeley [email protected]

, Rolf Turner [email protected]

and Ege Rubak [email protected]

See Also

dpois, gauss.hermite.

Examples

dmixpois(7, 10, 1, invlink = I)
  dpois(7, 10)

  pmixpois(7, log(10), 0.2)
  ppois(7, 10)

  qmixpois(0.95, log(10), 0.2)
  qpois(0.95, 10)

  x <- rmixpois(100, log(10), log(1.2))
  mean(x)
  var(x)

---

Extract the Domain of any Spatial Object

Description

Given a spatial object such as a point pattern, in any number of dimensions, this function extracts the spatial domain in which the object is defined.

Usage

## S3 method for class 'rmhmodel'
domain(X, ...)

Arguments

X	A spatial object such as a point pattern (in any number of dimensions), line segment pattern or pixel image.
...	Extra arguments. They are ignored by all the methods listed here.

Details

The function domain is generic.

For a spatial object X in any number of dimensions, domain(X) extracts the spatial domain in which X is defined.

For a two-dimensional object X, typically domain(X) is the same as Window(X).

Exceptions occur for methods related to linear networks.

Value

A spatial object representing the domain of X. Typically a window (object of class "owin"), a three-dimensional box ("box3"), a multidimensional box ("boxx") or a linear network ("linnet").

Author(s)

Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

See Also

domain, domain.quadratcount, domain.ppm, domain.quadrattest, domain.lpp. Window, Frame.

Examples

domain(rmhmodel(cif='poisson', par=list(beta=1), w=square(2)))

---

Pakes distribution

Description

Probability density, cumulative distribution function, quantile function, and random generation for the Pakes distribution.

Usage

dpakes(x, zeta)
ppakes(q, zeta)
qpakes(p, zeta)
rpakes(n, zeta)

Arguments

x, q	Numeric vector of quantiles.
p	Numeric vector of probabilities
n	Number of observations.
zeta	Mean of distribution. A single, non-negative, numeric value.

Details

These functions concern the probability distribution of the random variable

$X = \sum_{n=1}^\infty \prod_{j=1}^n U_j^{1/\zeta}$

where $U_1, U_2, \ldots$ are independent random variables uniformly distributed on $[0,1]$ and $\zeta$ is a parameter.

This distribution arises in many contexts. For example, for a homogeneous Poisson point process in two-dimensional space with intensity $\lambda$, the standard Gaussian kernel estimator of intensity with bandwidth $\sigma$, evaluated at any fixed location $u$, has the same distribution as $(\lambda/\zeta) X$ where $\zeta = 2 \pi \lambda\sigma^2$.

Following the usual convention, dpakes computes the probability density, ppakes the cumulative distribution function, and qpakes the quantile function, and rpakes generates random variates with this distribution.

The computation is based on a recursive integral equation for the cumulative distribution function, due to Professor Tony Pakes, presented in Baddeley, Moller and Pakes (2008). The solution uses the fact that the random variable satisfies the distributional equivalence

$X \equiv U^{1/\zeta} (1 + X)$

where $U$ is uniformly distributed on $[0,1]$ and independent of $X$.

Value

A numeric vector.

Author(s)

Adrian Baddeley.

References

Baddeley, A., Moller, J. and Pakes, A.G. (2008) Properties of residuals for spatial point processes, Annals of the Institute of Statistical Mathematics 60, 627–649.

Examples

curve(dpakes(x, 1.5), to=4)
 rpakes(3, 1.5)

---

Apply Expansion Rule

Description

Applies an expansion rule to a window.

Usage

expand.owin(W, ...)

Arguments

W	A window.
...	Arguments passed to rmhexpand to determine an expansion rule.

Details

The argument W should be a window (an object of class "owin").

This command applies the expansion rule specified by the arguments ... to the window W, yielding another window.

The arguments ... are passed to rmhexpand to determine the expansion rule.

For other transformations of the scale, location and orientation of a window, see shift, affine and rotate.

Value

A window (object of class "owin").

Author(s)

Adrian Baddeley [email protected]

and Rolf Turner [email protected]

See Also

rmhexpand about expansion rules.

shift, rotate, affine for other types of manipulation.

Examples

expand.owin(square(1), 9)
   expand.owin(square(1), distance=0.5)
   expand.owin(letterR, length=2)
   expand.owin(letterR, distance=0.1)

---

Gauss-Hermite Quadrature Approximation to Expectation for Normal Distribution

Description

Calculates an approximation to the expected value of any function of a normally-distributed random variable, using Gauss-Hermite quadrature.

Usage

gauss.hermite(f, mu = 0, sd = 1, ..., order = 5)

Arguments

f	The function whose moment should be approximated.
mu	Mean of the normal distribution.
sd	Standard deviation of the normal distribution.
...	Additional arguments passed to f.
order	Number of quadrature points in the Gauss-Hermite quadrature approximation. A small positive integer.

Details

This algorithm calculates the approximate expected value of f(Z) when Z is a normally-distributed random variable with mean mu and standard deviation sd. The expected value is an integral with respect to the Gaussian density; this integral is approximated using Gauss-Hermite quadrature.

The argument f should be a function in the R language whose first argument is the variable Z. Additional arguments may be passed through .... The value returned by f may be a single numeric value, a vector, or a matrix. The values returned by f for different values of Z must have compatible dimensions.

The result is a weighted average of several values of f.

Value

Numeric value, vector or matrix.

Author(s)

Adrian Baddeley [email protected]

, Rolf Turner [email protected]

and Ege Rubak [email protected].

Examples

gauss.hermite(function(x) x^2, 3, 1)

---

Recognise Stationary and Poisson Point Process Models

Description

Given a point process model (either a model that has been fitted to data, or a model specified by its parameters), determine whether the model is a stationary point process, and whether it is a Poisson point process.

Usage

is.stationary(x)
## S3 method for class 'rmhmodel'
is.stationary(x)

is.poisson(x)
## S3 method for class 'rmhmodel'
is.poisson(x)

Arguments

x	A fitted spatial point process model (object of class "ppm", "kppm", "lppm", "dppm" or "slrm") or a specification of a Gibbs point process model (object of class "rmhmodel") or a similar object.

Details

The argument x represents a fitted spatial point process model or a similar object.

is.stationary(x) returns TRUE if x represents a stationary point process, and FALSE if not.

is.poisson(x) returns TRUE if x represents a Poisson point process, and FALSE if not.

The functions is.stationary and is.poisson are generic, with methods for the classes "ppm" (Gibbs point process models), "kppm" (cluster or Cox point process models), "slrm" (spatial logistic regression models) and "rmhmodel" (model specifications for the Metropolis-Hastings algorithm). Additionally is.stationary has a method for classes "detpointprocfamily" and "dppm" (both determinantal point processes) and is.poisson has a method for class "interact" (interaction structures for Gibbs models).

is.poisson.kppm will return FALSE, unless the model x is degenerate: either x has zero intensity so that its realisations are empty with probability 1, or it is a log-Gaussian Cox process where the log intensity has zero variance.

is.poisson.slrm will always return TRUE, by convention.

Value

A logical value.

Author(s)

Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

See Also

is.marked to determine whether a model is a marked point process.

Examples

m <- rmhmodel(cif='strauss', par=list(beta=10, gamma=0.1, r=1))
  is.stationary(m)
  is.poisson(m)
  is.poisson(rmhmodel(cif='strauss', par=list(beta=10, gamma=1, r=1)))

---

Resample a Point Pattern by Resampling Quadrats

Description

Given a point pattern dataset, create a resampled point pattern by dividing the window into rectangular quadrats and randomly resampling the list of quadrats.

Usage

quadratresample(X, nx, ny=nx, ...,
                replace = FALSE, nsamples = 1,
                verbose = (nsamples > 1))

Arguments

X	A point pattern dataset (object of class "ppp").
nx, ny	Numbers of quadrats in the $x$ and $y$ directions.
...	Ignored.
replace	Logical value. Specifies whether quadrats should be sampled with or without replacement.
nsamples	Number of randomised point patterns to be generated.
verbose	Logical value indicating whether to print progress reports.

Details

This command implements a very simple bootstrap resampling procedure for spatial point patterns X.

The dataset X must be a point pattern (object of class "ppp") and its observation window must be a rectangle.

The window is first divided into N = nx * ny rectangular tiles (quadrats) of equal size and shape. To generate one resampled point pattern, a random sample of N quadrats is selected from the list of N quadrats, with replacement (if replace=TRUE) or without replacement (if replace=FALSE). The $i$th quadrat in the original dataset is then replaced by the $i$th sampled quadrat, after the latter is shifted so that it occupies the correct spatial position. The quadrats are then reconstituted into a point pattern inside the same window as X.

If replace=FALSE, this procedure effectively involves a random permutation of the quadrats. The resulting resampled point pattern has the same number of points as X. If replace=TRUE, the number of points in the resampled point pattern is random.

Value

A point pattern (if nsamples = 1) or a list of point patterns (if nsamples > 1).

Author(s)

Adrian Baddeley [email protected] and Rolf Turner [email protected]

See Also

quadrats, quadratcount.

See varblock to estimate the variance of a summary statistic by block resampling.

Examples

quadratresample(bei, 6, 3)

---

Alternating Gibbs Sampler for Multitype Point Processes

Description

Simulate a realisation of a point process model using the alternating Gibbs sampler.

Usage

rags(model, ..., ncycles = 100)

Arguments

model	Data specifying some kind of point process model.
...	Additional arguments passed to other code.
ncycles	Number of cycles of the alternating Gibbs sampler that should be performed.

Details

The Alternating Gibbs Sampler for a multitype point process is an iterative simulation procedure. Each step of the sampler updates the pattern of points of a particular type i, by drawing a realisation from the conditional distribution of points of type i given the points of all other types. Successive steps of the sampler update the points of type 1, then type 2, type 3, and so on.

This is an experimental implementation which currently works only for multitype hard core processes (see MultiHard) in which there is no interaction between points of the same type.

The argument model should be an object describing a point process model. At the moment, the only permitted format for model is of the form list(beta, hradii) where beta gives the first order trend and hradii is the matrix of interaction radii. See ragsMultiHard for full details.

Value

A point pattern (object of class "ppp").

Author(s)

Adrian Baddeley [email protected]

See Also

ragsMultiHard, ragsAreaInter

Examples

mo <- list(beta=c(30, 20),
             hradii = 0.05 * matrix(c(0,1,1,0), 2, 2))
  rags(mo, ncycles=10)

---

Alternating Gibbs Sampler for Area-Interaction Process

Description

Generate a realisation of the area-interaction process using the alternating Gibbs sampler. Applies only when the interaction parameter $eta$ is greater than 1.

Usage

ragsAreaInter(beta, eta, r, ...,
                   win = NULL, bmax = NULL, periodic = FALSE, ncycles = 100)

Arguments

beta	First order trend. A number, a pixel image (object of class "im"), or a function(x,y).
eta	Interaction parameter (canonical form) as described in the help for AreaInter. A number greater than 1.
r	Disc radius in the model. A number greater than 1.
...	Additional arguments for beta if it is a function.
win	Simulation window. An object of class "owin". (Ignored if beta is a pixel image.)
bmax	Optional. The maximum possible value of beta, or a number larger than this.
periodic	Logical value indicating whether to treat opposite sides of the simulation window as being the same, so that points close to one side may interact with points close to the opposite side. Feasible only when the window is a rectangle.
ncycles	Number of cycles of the alternating Gibbs sampler to be performed.

Details

This function generates a simulated realisation of the area-interaction process (see AreaInter) using the alternating Gibbs sampler (see rags).

It exploits a mathematical relationship between the (unmarked) area-interaction process and the two-type hard core process (Baddeley and Van Lieshout, 1995; Widom and Rowlinson, 1970). This relationship only holds when the interaction parameter eta is greater than 1 so that the area-interaction process is clustered.

The parameters beta,eta are the canonical parameters described in the help for AreaInter. The first order trend beta may be a constant, a function, or a pixel image.

The simulation window is determined by beta if it is a pixel image, and otherwise by the argument win (the default is the unit square).

Value

A point pattern (object of class "ppp").

Author(s)

Adrian Baddeley [email protected].

References

Baddeley, A.J. and Van Lieshout, M.N.M. (1995). Area-interaction point processes. Annals of the Institute of Statistical Mathematics 47 (1995) 601–619.

Widom, B. and Rowlinson, J.S. (1970). New model for the study of liquid-vapor phase transitions. The Journal of Chemical Physics 52 (1970) 1670–1684.

See Also

rags, ragsMultiHard

AreaInter

Examples

plot(ragsAreaInter(100, 2, 0.07, ncycles=15))

---

Alternating Gibbs Sampler for Multitype Hard Core Process

Description

Generate a realisation of the multitype hard core point process using the alternating Gibbs sampler.

Usage

ragsMultiHard(beta, hradii, ..., types=NULL, bmax = NULL,
              periodic=FALSE, ncycles = 100)

Arguments

beta	First order trend. A numeric vector, a pixel image, a function, a list of functions, or a list of pixel images.
hradii	Matrix of hard core radii between each pair of types. Diagonal entries should be 0 or NA.
types	Vector of all possible types for the multitype point pattern.
...	Arguments passed to rmpoispp when generating random points.
bmax	Optional upper bound on beta.
periodic	Logical value indicating whether to measure distances in the periodic sense, so that opposite sides of the (rectangular) window are treated as identical.
ncycles	Number of cycles of the sampler to be performed.

Details

The Alternating Gibbs Sampler for a multitype point process is an iterative simulation procedure. Each step of the sampler updates the pattern of points of a particular type i, by drawing a realisation from the conditional distribution of points of type i given the points of all other types. Successive steps of the sampler update the points of type 1, then type 2, type 3, and so on.

This is an experimental implementation which currently works only for multitype hard core processes (see MultiHard) in which there is no interaction between points of the same type, and for the area-interaction process (see ragsAreaInter).

The argument beta gives the first order trend for each possible type of point. It may be a single number, a numeric vector, a function(x,y), a pixel image, a list of functions, a function(x,y,m), or a list of pixel images.

The argument hradii is the matrix of hard core radii between each pair of possible types of points. Two points of types i and j respectively are forbidden to lie closer than a distance hradii[i,j] apart. The diagonal of this matrix must contain NA or 0 values, indicating that there is no hard core constraint applying between points of the same type.

Value

A point pattern (object of class "ppp").

Author(s)

Adrian Baddeley [email protected]

See Also

rags, ragsAreaInter

Examples

b <- c(30,20)
  h <- 0.05 * matrix(c(0,1,1,0), 2, 2)
  ragsMultiHard(b, h, ncycles=10)
  ragsMultiHard(b, h, ncycles=5, periodic=TRUE)

---

Simulate Neyman-Scott Point Process with Cauchy cluster kernel

Description

Generate a random point pattern, a simulated realisation of the Neyman-Scott process with Cauchy cluster kernel.

Usage

rCauchy(kappa, scale, mu, win = square(1),
            nsim=1, drop=TRUE,
            ...,
            algorithm=c("BKBC", "naive"),
            nonempty=TRUE, 
            thresh = 0.001, poisthresh=1e-6,
            expand = NULL,
            saveparents=FALSE, saveLambda=FALSE,
            kappamax=NULL, mumax=NULL)

Arguments

kappa	Intensity of the Poisson process of cluster centres. A single positive number, a function, or a pixel image.
scale	Scale parameter for cluster kernel. Determines the size of clusters. A single positive number, in the same units as the spatial coordinates.
mu	Mean number of points per cluster (a single positive number) or reference intensity for the cluster points (a function or a pixel image).
win	Window in which to simulate the pattern. An object of class "owin" or something acceptable to as.owin.
nsim	Number of simulated realisations to be generated.
drop	Logical. If nsim=1 and drop=TRUE (the default), the result will be a point pattern, rather than a list containing a point pattern.
...	Passed to clusterfield to control the image resolution when saveLambda=TRUE, and to clusterradius when expand is missing or NULL.
algorithm	String (partially matched) specifying the simulation algorithm. See Details.
nonempty	Logical. If TRUE (the default), a more efficient algorithm is used, in which parents are generated conditionally on having at least one offspring point. If FALSE, parents are generated even if they have no offspring. Both choices are valid; the default is recommended unless you need to simulate all the parent points for some other purpose.
thresh	Threshold relative to the cluster kernel value at the origin (parent location) determining when the cluster kernel will be treated as zero for simulation purposes. Will be overridden by argument expand if that is given.
poisthresh	Numerical threshold below which the model will be treated as a Poisson process. See Details.
expand	Window expansion distance. A single number. The distance by which the original window will be expanded in order to generate parent points. Has a sensible default, determined by calling clusterradius with the numeric threshold value given in thresh.
saveparents	Logical value indicating whether to save the locations of the parent points as an attribute.
saveLambda	Logical. If TRUE then the random intensity corresponding to the simulated parent points will also be calculated and saved, and returns as an attribute of the point pattern.
kappamax	Optional. Numerical value which is an upper bound for the values of kappa, when kappa is a pixel image or a function.
mumax	Optional. Numerical value which is an upper bound for the values of mu, when mu is a pixel image or a function.

Details

This algorithm generates a realisation of the Neyman-Scott process with Cauchy cluster kernel, inside the window win.

The process is constructed by first generating a Poisson point process of “parent” points with intensity kappa. Then each parent point is replaced by a random cluster of points, the number of points in each cluster being random with a Poisson (mu) distribution, and the points being placed independently and uniformly according to a Cauchy kernel.

Note that, for correct simulation of the model, the parent points are not restricted to lie inside the window win; the parent process is effectively the uniform Poisson process on the infinite plane.

The algorithm can also generate spatially inhomogeneous versions of the cluster process:

• The parent points can be spatially inhomogeneous. If the argument kappa is a function(x,y) or a pixel image (object of class "im"), then it is taken as specifying the intensity function of an inhomogeneous Poisson process that generates the parent points.

• The offspring points can be inhomogeneous. If the argument mu is a function(x,y) or a pixel image (object of class "im"), then it is interpreted as the reference density for offspring points, in the sense of Waagepetersen (2006).

When the parents are homogeneous (kappa is a single number) and the offspring are inhomogeneous (mu is a function or pixel image), the model can be fitted to data using kppm.

If the pair correlation function of the model is very close to that of a Poisson process, deviating by less than poisthresh, then the model is approximately a Poisson process, and will be simulated as a Poisson process with intensity kappa * mu, using rpoispp. This avoids computations that would otherwise require huge amounts of memory.

Value

A point pattern (an object of class "ppp") if nsim=1, or a list of point patterns if nsim > 1.

Additionally, some intermediate results of the simulation are returned as attributes of this point pattern (see rNeymanScott). Furthermore, the simulated intensity function is returned as an attribute "Lambda", if saveLambda=TRUE.

Simulation Algorithm

Two simulation algorithms are implemented.

• The naive algorithm generates the cluster process by directly following the description given above. First the window win is expanded by a distance equal to expand. Then the parent points are generated in the expanded window according to a Poisson process with intensity kappa. Then each parent point is replaced by a finite cluster of offspring points as described above. The naive algorithm is used if algorithm="naive" or if nonempty=FALSE.

• The BKBC algorithm, proposed by Baddeley and Chang (2023), is a modification of the algorithm of Brix and Kendall (2002). Parents are generated in the infinite plane, subject to the condition that they have at least one offspring point inside the window win. The BKBC algorithm is used when algorithm="BKBC" (the default) and nonempty=TRUE (the default).

The naive algorithm becomes very slow when scale is large, while the BKBC algorithm is uniformly fast (Baddeley and Chang, 2023).

If saveparents=TRUE, then the simulated point pattern will have an attribute "parents" containing the coordinates of the parent points, and an attribute "parentid" mapping each offspring point to its parent.

If nonempty=TRUE (the default), then parents are generated subject to the condition that they have at least one offspring point in the window win. nonempty=FALSE, then parents without offspring will be included; this option is not available in the BKBC algorithm.

Note that if kappa is a pixel image, its domain must be larger than the window win. This is because an offspring point inside win could have its parent point lying outside win. In order to allow this, the naive simulation algorithm first expands the original window win by a distance equal to expand and generates the Poisson process of parent points on this larger window. If kappa is a pixel image, its domain must contain this larger window.

If the pair correlation function of the model is very close to that of a Poisson process, with maximum deviation less than poisthresh, then the model is approximately a Poisson process. This is detected by the naive algorithm which then simulates a Poisson process with intensity kappa * mu, using rpoispp. This avoids computations that would otherwise require huge amounts of memory.

Fitting cluster models to data

The Cauchy cluster model with homogeneous parents (i.e. where kappa is a single number) where the offspring are either homogeneous or inhomogeneous (mu is a single number, a function or pixel image) can be fitted to point pattern data using kppm, or fitted to the inhomogeneous $K$ function using cauchy.estK or cauchy.estpcf.

Currently spatstat does not support fitting the Cauchy cluster process model with inhomogeneous parents.

A Cauchy cluster process model fitted by kppm can be simulated automatically using simulate.kppm (which invokes rCauchy to perform the simulation).

Author(s)

Original algorithm by Abdollah Jalilian and Rasmus Waagepetersen. Adapted for spatstat by Adrian Baddeley [email protected]. Brix-Kendall-Baddeley-Chang algorithm implemented by Adrian Baddeley [email protected] and Ya-Mei Chang [email protected].

References

Baddeley, A. and Chang, Y.-M. (2023) Robust algorithms for simulating cluster point processes. Journal of Statistical Computation and Simulation. In Press. DOI 10.1080/00949655.2023.2166045.

Brix, A. and Kendall, W.S. (2002) Simulation of cluster point processes without edge effects. Advances in Applied Probability 34, 267–280.

Ghorbani, M. (2013) Cauchy cluster process. Metrika 76, 697-706.

Jalilian, A., Guan, Y. and Waagepetersen, R. (2013) Decomposition of variance for spatial Cox processes. Scandinavian Journal of Statistics 40, 119-137.

Waagepetersen, R. (2007) An estimating function approach to inference for inhomogeneous Neyman-Scott processes. Biometrics 63, 252–258.

See Also

rpoispp, rMatClust, rThomas, rVarGamma, rNeymanScott, rGaussPoisson.

For fitting the model, see kppm, clusterfit.

Examples

# homogeneous
 X <- rCauchy(30, 0.01, 5)
 # inhomogeneous
 ff <- function(x,y){ exp(2 - 3 * abs(x)) }
 Z <- as.im(ff, W= owin())
 Y <- rCauchy(50, 0.01, Z)
 YY <- rCauchy(ff, 0.01, 5)

---

Simulate Baddeley-Silverman Cell Process

Description

Generates a random point pattern, a simulated realisation of the Baddeley-Silverman cell process model.

Usage

rcell(win=square(1), nx=NULL, ny=nx, ..., dx=NULL, dy=dx,
       N=10, nsim=1, drop=TRUE)

Arguments

win	A window. An object of class owin, or data in any format acceptable to as.owin().
nx	Number of columns of cells in the window. Incompatible with dx.
ny	Number of rows of cells in the window. Incompatible with dy.
...	Ignored.
dx	Width of the cells. Incompatible with nx.
dy	Height of the cells. Incompatible with ny.
N	Integer. Distributional parameter: the maximum number of random points in each cell. Passed to rcellnumber.
nsim	Number of simulated realisations to be generated.
drop	Logical. If nsim=1 and drop=TRUE (the default), the result will be a point pattern, rather than a list containing a point pattern.

Details

This function generates a simulated realisation of the “cell process” (Baddeley and Silverman, 1984), a random point process with the same second-order properties as the uniform Poisson process. In particular, the $K$ function of this process is identical to the $K$ function of the uniform Poisson process (aka Complete Spatial Randomness). The same holds for the pair correlation function and all other second-order properties. The cell process is a counterexample to the claim that the $K$ function completely characterises a point pattern.

A cell process is generated by dividing space into equal rectangular tiles. In each tile, a random number of random points is placed. By default, there are either $0$, $1$ or $10$ points, with probabilities $1/10$, $8/9$ and $1/90$ respectively. The points within a tile are independent and uniformly distributed in that tile, and the numbers of points in different tiles are independent random integers.

The tile width is determined either by the number of columns nx or by the horizontal spacing dx. The tile height is determined either by the number of rows ny or by the vertical spacing dy. The cell process is then generated in these tiles. The random numbers of points are generated by rcellnumber.

Some of the resulting random points may lie outside the window win: if they do, they are deleted. The result is a point pattern inside the window win.

Value

A point pattern (an object of class "ppp") if nsim=1, or a list of point patterns if nsim > 1.

Author(s)

Adrian Baddeley [email protected] and Rolf Turner [email protected]

References

Baddeley, A.J. and Silverman, B.W. (1984) A cautionary example on the use of second-order methods for analyzing point patterns. Biometrics 40, 1089-1094.

See Also

rcellnumber, rstrat, rsyst, runifpoint, Kest

Examples

X <- rcell(nx=15)
  plot(X)
  if(require(spatstat.explore)) {
     plot(Kest(X))
  }

---

Generate Random Numbers of Points for Cell Process

Description

Generates random integers for the Baddeley-Silverman counterexample.

Usage

rcellnumber(n, N = 10, mu=1)

Arguments

n	Number of random integers to be generated.
N	Distributional parameter: the largest possible value (when mu <= 1). An integer greater than 1.
mu	Mean of the distribution (equals the variance). Any positive real number.

Details

If mu = 1 (the default), this function generates random integers which have mean and variance equal to 1, but which do not have a Poisson distribution. The random integers take the values $0$, $1$ and $N$ with probabilities $1/N$, $(N-2)/(N-1)$ and $1/(N(N-1))$ respectively. See Baddeley and Silverman (1984).

If mu is another positive number, the random integers will have mean and variance equal to mu. They are obtained by generating the one-dimensional counterpart of the cell process and counting the number of points in the interval from 0 to mu. The maximum possible value of each random integer is N * ceiling(mu).

Value

An integer vector of length n.

Author(s)

Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

References

Baddeley, A.J. and Silverman, B.W. (1984) A cautionary example on the use of second-order methods for analyzing point patterns. Biometrics 40, 1089-1094.

See Also

rcell

Examples

rcellnumber(30, 3)

---

Simulate Cluster Process using Brix-Kendall Algorithm or Modifications

Description

Generates simulated realisations of a stationary Neyman-Scott cluster point process, using the Brix-Kendall (2002) algorithm or various modifications proposed by Baddeley and Chang (2023). For advanced research use.

Usage

rclusterBKBC(clusters="Thomas",
   kappa, mu, scale,
   ...,
   W = unit.square(),
   nsim = 1, drop = TRUE,
   best = FALSE,
   external = c("BK", "superBK", "border"),
   internal = c("dominating", "naive"),
   inflate = 1,
   psmall = 1e-04,
   use.inverse=TRUE,
   use.special=TRUE,
   integralmethod=c("quadrature", "trapezoid"),
   verbose = TRUE, warn=TRUE)

Arguments

clusters	Character string (partially matched) specifying the cluster process. Current options include "Thomas", "MatClust", "Cauchy" and "VarGamma".
kappa	Intensity of the parent process. A nonnegative number.
mu	Mean number of offspring per parent. A nonnegative number.
scale	Cluster scale. Interpretation depends on the model.
...	Additional arguments controlling the shape of the cluster kernel, if any.
W	Window in which the simulation should be generated. An object of class "owin".
nsim	The number of simulated point patterns to be generated. A positive integer.
drop	Logical value. If nsim=1 and drop=TRUE (the default), the result will be a point pattern, rather than a list containing a point pattern.
best	Logical value. If best=TRUE, the code will choose the fastest algorithm. If best=FALSE (the default), the algorithm will be specified by the other arguments external and internal. See Details.
external	Algorithm to be used to generate parent points which lie outside the bounding window. See Details.
internal	Algorithm to be used to generate parent points which lie inside the bounding window. See Details.
inflate	Numerical value determining the position of the bounding window. See Details.
psmall	Threshold of small probability for use in the algorithm.
use.inverse	Logical value specifying whether to compute the inverse function analytically, if possible (use.inverse=TRUE, the default) or by numerical root-finding (use.inverse=FALSE). This is mainly for checking validity of code.
use.special	Logical value specifying whether to use efficient special code (if available) to generate the simulations (use.special=TRUE, the default) or to use generic code (use.special=FALSE). This is mainly for checking validity of code.
integralmethod	Character string (partially matched) specifying how to perform numerical computation of integrals when required. This argument is passed to indefinteg. The default integralmethod="quadrature" is accurate but can be slow. Faster, but possibly less accurate, integration can be performed by setting integralmethod="trapezoid".
verbose	Logical value specifying whether to print detailed information about the simulation algorithm during execution.
warn	Logical value specifying whether to issue a warning if the number of random proposal points is very large.

Details

This function is intended for advanced research use. It implements the algorithm of Brix and Kendall (2002) for generating simulated realisations of a stationary Neyman-Scott process, and various modifications of this algorithm proposed in Baddeley and Chang (2023). It is an alternative to rNeymanScott.

The function supports the following models:

• clusters="Thomas": the (modified) Thomas cluster process which can also be simulated by rThomas.

• clusters="MatClust": the Matern cluster process which can also be simulated by rMatClust.

• clusters="Cauchy": the Cauchy cluster process which can also be simulated by rCauchy.

• clusters="VarGamma": the variance-gamma cluster process which can also be simulated by rVarGamma.

• any other Poisson cluster process models that may be recognised by kppm.

By default, the code executes the original Brix-Kendall algorithm described in Sections 2.3 and 3.1 of Brix and Kendall (2002).

Modifications of this algorithm, proposed in Baddeley and Chang (2023), can be selected using the arguments external and internal, or best.

If best=TRUE, the code will choose the algorithm that would run fastest with the given parameters. If best=FALSE (the default), the choice of algorithm is determined by the arguments external and internal.

First the window W is enclosed in a disc D and Monte Carlo proposal densities are defined with reference to D as described in Brix and Kendall (2002). Then D is inflated by the scale factor inflate to produce a larger disc E (by default inflate=1 implying E=D). Then the parent points of the clusters are generated, possibly using different mechanisms inside and outside E.

The argument external determines the algorithm for generating parent points outside E.

• If external="BK" (the default), proposed parents outside E will be generated from a dominating point process as described in Section 3.1 of Brix and Kendall (2002). These points will be thinned to obtain the correct intensity of parent points. For each accepted parent, offspring points are generated inside D, subject to the condition that the parent has at least one offspring inside D. Offspring points are subsequently clipped to the true window W.

• If external="superBK", proposed parents will initially be generated from a process that dominates the dominating point process as described in Baddeley and Chang (2023). These proposals will then be thinned to obtain the correct intensity of the dominating process, then thinned again to obtain the correct intensity of parent points. This procedure reduces computation time when scale is large. For each accepted parent, offspring points are generated inside D, subject to the condition that the parent has at least one offspring inside D. Offspring points are subsequently clipped to the true window W.

• If external="border" then proposed parents will be generated with uniform intensity in a border region surrounding the disc D. For each proposed parent, offspring points are generated in the entire plane according to the cluster offspring distribution, without any restriction. Offspring points are subsequently clipped to the true window W. This is the technique currently used in rNeymanScott.

The argument internal determines the algorithm for generating proposed parent points inside E.

• If internal="dominating", parent points in E are generated according to the dominating point process described in Sections 2.3 and 3.1 of Brix and Kendall (2002), and then thinned to obtain the correct intensity of parent points. For each accepted parent, offspring points are generated inside D, subject to the condition that the parent has at least one offspring inside D. Offspring points are subsequently clipped to the true window W.

• If internal="naive", parent points in E are generated with uniform intensity inside E and are not thinned. For each proposed parent, offspring points are generated in the entire plane according to the cluster offspring distribution, without any restriction. Offspring points are subsequently clipped to the true window W. This is the technique currently used in rNeymanScott.

If warn=TRUE, then a warning will be issued if the number of random proposal points (proposed parents and proposed offspring) is very large. The threshold is spatstat.options("huge.npoints"). This warning has no consequences, but it helps to trap a number of common problems.

Value

A point pattern, or a list of point patterns.

If nsim=1 and drop=TRUE, the result is a point pattern (an object of class "ppp").

Otherwise, the result is a list of nsim point patterns, and also belongs to the class "solist".

Author(s)

Adrian Baddeley [email protected] and Ya-Mei Chang [email protected].

References

Baddeley, A. and Chang, Y.-M. (2023) Robust algorithms for simulating cluster point processes. Journal of Statistical Computation and Simulation. In Press. DOI 10.1080/00949655.2023.2166045.

Brix, A. and Kendall, W.S. (2002) Simulation of cluster point processes without edge effects. Advances in Applied Probability 34, 267–280.

See Also

rNeymanScott, rMatClust, rThomas, rCauchy, rVarGamma

Examples

Y <- rclusterBKBC("Thomas", 10,5,0.2)
  Y
  Z <- rclusterBKBC("VarGamma", 10,5,0.2,
          nu=-1/4,
          internal="naive", external="super",
          verbose=FALSE)

---

Perfect Simulation of the Diggle-Gates-Stibbard Process

Description

Generate a random pattern of points, a simulated realisation of the Diggle-Gates-Stibbard process, using a perfect simulation algorithm.

Usage

rDGS(beta, rho, W = owin(), expand=TRUE, nsim=1, drop=TRUE)

Arguments

beta	intensity parameter (a positive number).
rho	interaction range (a non-negative number).
W	window (object of class "owin") in which to generate the random pattern.
expand	Logical. If FALSE, simulation is performed in the window W, which must be rectangular. If TRUE (the default), simulation is performed on a larger window, and the result is clipped to the original window W. Alternatively expand can be an object of class "rmhexpand" (see rmhexpand) determining the expansion method.
nsim	Number of simulated realisations to be generated.
drop	Logical. If nsim=1 and drop=TRUE (the default), the result will be a point pattern, rather than a list containing a point pattern.

Details

This function generates a realisation of the Diggle-Gates-Stibbard point process in the window W using a ‘perfect simulation’ algorithm.

Diggle, Gates and Stibbard (1987) proposed a pairwise interaction point process in which each pair of points separated by a distance $d$ contributes a factor $e(d)$ to the probability density, where

$e(d) = \sin^2\left(\frac{\pi d}{2\rho}\right)$

for $d < \rho$, and $e(d)$ is equal to 1 for $d \ge \rho$.

The simulation algorithm used to generate the point pattern is ‘dominated coupling from the past’ as implemented by Berthelsen and Moller (2002, 2003). This is a ‘perfect simulation’ or ‘exact simulation’ algorithm, so called because the output of the algorithm is guaranteed to have the correct probability distribution exactly (unlike the Metropolis-Hastings algorithm used in rmh, whose output is only approximately correct).

There is a tiny chance that the algorithm will run out of space before it has terminated. If this occurs, an error message will be generated.

Value

If nsim = 1, a point pattern (object of class "ppp"). If nsim > 1, a list of point patterns.

Author(s)

Adrian Baddeley [email protected], based on original code for the Strauss process by Kasper Klitgaard Berthelsen.

References

Berthelsen, K.K. and Moller, J. (2002) A primer on perfect simulation for spatial point processes. Bulletin of the Brazilian Mathematical Society 33, 351-367.

Berthelsen, K.K. and Moller, J. (2003) Likelihood and non-parametric Bayesian MCMC inference for spatial point processes based on perfect simulation and path sampling. Scandinavian Journal of Statistics 30, 549-564.

Diggle, P.J., Gates, D.J., and Stibbard, A. (1987) A nonparametric estimator for pairwise-interaction point processes. Biometrika 74, 763 – 770. Scandinavian Journal of Statistics 21, 359–373.

Moller, J. and Waagepetersen, R. (2003). Statistical Inference and Simulation for Spatial Point Processes. Chapman and Hall/CRC.

See Also

rmh, DiggleGatesStibbard.

rStrauss, rHardcore, rStraussHard, rDiggleGratton, rPenttinen.

Examples

X <- rDGS(50, 0.05)
   Z <- rDGS(50, 0.03, nsim=2)

---

Perfect Simulation of the Diggle-Gratton Process

Description

Generate a random pattern of points, a simulated realisation of the Diggle-Gratton process, using a perfect simulation algorithm.

Usage

rDiggleGratton(beta, delta, rho, kappa=1, W = owin(),
                 expand=TRUE, nsim=1, drop=TRUE)

Arguments

beta	intensity parameter (a positive number).
delta	hard core distance (a non-negative number).
rho	interaction range (a number greater than delta).
kappa	interaction exponent (a non-negative number).
W	window (object of class "owin") in which to generate the random pattern. Currently this must be a rectangular window.
expand	Logical. If FALSE, simulation is performed in the window W, which must be rectangular. If TRUE (the default), simulation is performed on a larger window, and the result is clipped to the original window W. Alternatively expand can be an object of class "rmhexpand" (see rmhexpand) determining the expansion method.
nsim	Number of simulated realisations to be generated.
drop	Logical. If nsim=1 and drop=TRUE (the default), the result will be a point pattern, rather than a list containing a point pattern.

Details

This function generates a realisation of the Diggle-Gratton point process in the window W using a ‘perfect simulation’ algorithm.

Diggle and Gratton (1984, pages 208-210) introduced the pairwise interaction point process with pair potential $h(t)$ of the form

$h(t) = \left( \frac{t-\delta}{\rho-\delta} \right)^\kappa \quad\quad \mbox{ if } \delta \le t \le \rho$

with $h(t) = 0$ for $t < \delta$ and $h(t) = 1$ for $t > \rho$. Here $\delta$, $\rho$ and $\kappa$ are parameters.

Note that we use the symbol $\kappa$ where Diggle and Gratton (1984) use $\beta$, since in spatstat we reserve the symbol $\beta$ for an intensity parameter.

The parameters must all be nonnegative, and must satisfy $\delta \le \rho$.

The simulation algorithm used to generate the point pattern is ‘dominated coupling from the past’ as implemented by Berthelsen and Moller (2002, 2003). This is a ‘perfect simulation’ or ‘exact simulation’ algorithm, so called because the output of the algorithm is guaranteed to have the correct probability distribution exactly (unlike the Metropolis-Hastings algorithm used in rmh, whose output is only approximately correct).

There is a tiny chance that the algorithm will run out of space before it has terminated. If this occurs, an error message will be generated.

Value

If nsim = 1, a point pattern (object of class "ppp"). If nsim > 1, a list of point patterns.

Author(s)

Adrian Baddeley [email protected]

based on original code for the Strauss process by Kasper Klitgaard Berthelsen.

References

Berthelsen, K.K. and Moller, J. (2002) A primer on perfect simulation for spatial point processes. Bulletin of the Brazilian Mathematical Society 33, 351-367.

Berthelsen, K.K. and Moller, J. (2003) Likelihood and non-parametric Bayesian MCMC inference for spatial point processes based on perfect simulation and path sampling. Scandinavian Journal of Statistics 30, 549-564.

Diggle, P.J. and Gratton, R.J. (1984) Monte Carlo methods of inference for implicit statistical models. Journal of the Royal Statistical Society, series B 46, 193 – 212.

Moller, J. and Waagepetersen, R. (2003). Statistical Inference and Simulation for Spatial Point Processes. Chapman and Hall/CRC.

See Also

rmh, rStrauss, rHardcore, rStraussHard, rDGS, rPenttinen.

For fitting the model, see DiggleGratton.

Examples

X <- rDiggleGratton(50, 0.02, 0.07)
   Z <- rDiggleGratton(50, 0.02, 0.07, 2, nsim=2)

---

Interaction Distance of a Point Process Model

Description

Computes the interaction distance of a point process model.

Usage

reach(x, ...)

  ## S3 method for class 'rmhmodel'
reach(x, ...)

Arguments

x	Either a fitted point process model (object of class "ppm"), an interpoint interaction (object of class "interact"), a fitted interpoint interaction (object of class "fii") or a point process model for simulation (object of class "rmhmodel").
...	Other arguments are ignored.

Details

The function reach computes the ‘interaction distance’ or ‘interaction range’ of a point process model.

The definition of the interaction distance depends on the type of point process model. This help page explains the interaction distance for a Gibbs point process. For other kinds of models, see reach.kppm and reach.dppm.

For a Gibbs point process model, the interaction distance is the shortest distance $D$ such that any two points in the process which are separated by a distance greater than $D$ do not interact with each other.

For example, the interaction range of a Strauss process (see Strauss or rStrauss) with parameters $\beta,\gamma,r$ is equal to $r$, unless $\gamma=1$ in which case the model is Poisson and the interaction range is $0$. The interaction range of a Poisson process is zero. The interaction range of the Ord threshold process (see OrdThresh) is infinite, since two points may interact at any distance apart.

The function reach is generic, with methods for the case where x is

• a fitted point process model (object of class "ppm", usually obtained from the model-fitting function ppm);

• an interpoint interaction structure (object of class "interact")

• a fitted interpoint interaction (object of class "fii")

• a point process model for simulation (object of class "rmhmodel"), usually obtained from rmhmodel.

Value

The interaction distance, or NA if this cannot be computed from the information given.

Other types of models

Methods for reach are also defined for point process models of class "kppm" and "dppm". Their technical definition is different from this one. See reach.kppm and reach.dppm.

Author(s)

Adrian Baddeley [email protected] and Rolf Turner [email protected]

See Also

reach.ppm

rmhmodel

See reach.kppm and reach.dppm for other types of point process models.

Examples

reach(rmhmodel(cif='poisson', par=list(beta=100)))
    # returns 0

    reach(rmhmodel(cif='strauss', par=list(beta=100, gamma=0.1, r=7)))
    # returns 7

    reach(rmhmodel(cif='sftcr', par=list(beta=100, sigma=1, kappa=0.7)))
    # returns Inf
    
    reach(rmhmodel(cif='multihard',
                   par=list(beta=c(10,10), hradii=matrix(c(1,3,3,1),2,2))))
    # returns 3

---

First Reciprocal Moment of the Truncated Poisson Distribution

Description

Computes the first reciprocal moment (first negative moment) of the truncated Poisson distribution (the Poisson distribution conditioned to have a nonzero value).

Usage

recipEnzpois(mu, exact=TRUE)

Arguments

mu	The mean of the original Poisson distribution. A single positive numeric value, or a vector of positive numbers.
exact	Logical value specifying whether to use the exact analytic formula if possible.

Details

This function calculates the expected value of $1/N$ given $N > 0$, where $N$ is a Poisson random variable with mean $\mu$.

If the library gsl is loaded, and if exact=TRUE (the default), then the calculation uses the exact analytic formula

$\nu = \frac{e^{-\mu}}{1- e^{-\mu}} \left( Ei(\mu) - \log \mu - \gamma \right)$

(see e.g. Grab and Savage, 1954) where $\nu$ is the desired reciprocal moment, and

$Ei(x) = \int_{-\infty}^x t e^{-t} dt$

is the first exponential integral, and $\gamma \approx 0.577$ is the Euler-Mascheroni constant.

If gsl is not loaded, or if exact=FALSE is specified, the value is computed approximately (and more slowly) by summing over the possible values of $N$ up to a finite limit.

Value

A single numerical value or a numeric vector.

Author(s)

Adrian Baddeley [email protected].

References

Grab, E.L. and Savage, I.R. (1954) Tables of the expected value of 1/X for positive Bernoulli and Poisson variables. Journal of the American Statistical Association 49, 169–177.

See Also

rpoisnonzero

Examples

if(require(gsl)) {
    v <- recipEnzpois(10)
    print(v)
  }
  recipEnzpois(10, exact=FALSE)

---

Simulate Gauss-Poisson Process

Description

Generate a random point pattern, a simulated realisation of the Gauss-Poisson Process.

Usage

rGaussPoisson(kappa, r, p2, win = owin(c(0,1),c(0,1)),
               ..., nsim=1, drop=TRUE)

Arguments

kappa	Intensity of the Poisson process of cluster centres. A single positive number, a function, or a pixel image.
r	Diameter of each cluster that consists of exactly 2 points.
p2	Probability that a cluster contains exactly 2 points.
win	Window in which to simulate the pattern. An object of class "owin" or something acceptable to as.owin.
...	Ignored.
nsim	Number of simulated realisations to be generated.
drop	Logical. If nsim=1 and drop=TRUE (the default), the result will be a point pattern, rather than a list containing a point pattern.

Details

This algorithm generates a realisation of the Gauss-Poisson point process inside the window win. The process is constructed by first generating a Poisson point process of parent points with intensity kappa. Then each parent point is either retained (with probability 1 - p2) or replaced by a pair of points at a fixed distance r apart (with probability p2). In the case of clusters of 2 points, the line joining the two points has uniform random orientation.

In this implementation, parent points are not restricted to lie in the window; the parent process is effectively the uniform Poisson process on the infinite plane.

Value

A point pattern (an object of class "ppp") if nsim=1, or a list of point patterns if nsim > 1.

Additionally, some intermediate results of the simulation are returned as attributes of the point pattern. See rNeymanScott.

Author(s)

Adrian Baddeley [email protected] and Rolf Turner [email protected]

See Also

rpoispp, rThomas, rMatClust, rNeymanScott

Examples

pp <- rGaussPoisson(30, 0.07, 0.5)

---

Perfect Simulation of the Hardcore Process

Description

Generate a random pattern of points, a simulated realisation of the Hardcore process, using a perfect simulation algorithm.

Usage

rHardcore(beta, R = 0, W = owin(), expand=TRUE, nsim=1, drop=TRUE)

Arguments

beta	intensity parameter (a positive number).
R	hard core distance (a non-negative number).
W	window (object of class "owin") in which to generate the random pattern. Currently this must be a rectangular window.
expand	Logical. If FALSE, simulation is performed in the window W, which must be rectangular. If TRUE (the default), simulation is performed on a larger window, and the result is clipped to the original window W. Alternatively expand can be an object of class "rmhexpand" (see rmhexpand) determining the expansion method.
nsim	Number of simulated realisations to be generated.
drop	Logical. If nsim=1 and drop=TRUE (the default), the result will be a point pattern, rather than a list containing a point pattern.

Details

This function generates a realisation of the Hardcore point process in the window W using a ‘perfect simulation’ algorithm.

The Hardcore process is a model for strong spatial inhibition. Two points of the process are forbidden to lie closer than R units apart. The Hardcore process is the special case of the Strauss process (see rStrauss) with interaction parameter $\gamma$ equal to zero.

The simulation algorithm used to generate the point pattern is ‘dominated coupling from the past’ as implemented by Berthelsen and Moller (2002, 2003). This is a ‘perfect simulation’ or ‘exact simulation’ algorithm, so called because the output of the algorithm is guaranteed to have the correct probability distribution exactly (unlike the Metropolis-Hastings algorithm used in rmh, whose output is only approximately correct).

There is a tiny chance that the algorithm will run out of space before it has terminated. If this occurs, an error message will be generated.

Value

If nsim = 1, a point pattern (object of class "ppp"). If nsim > 1, a list of point patterns.

Author(s)

Adrian Baddeley [email protected], based on original code for the Strauss process by Kasper Klitgaard Berthelsen.

References

Berthelsen, K.K. and Moller, J. (2002) A primer on perfect simulation for spatial point processes. Bulletin of the Brazilian Mathematical Society 33, 351-367.

Berthelsen, K.K. and Moller, J. (2003) Likelihood and non-parametric Bayesian MCMC inference for spatial point processes based on perfect simulation and path sampling. Scandinavian Journal of Statistics 30, 549-564.

Moller, J. and Waagepetersen, R. (2003). Statistical Inference and Simulation for Spatial Point Processes. Chapman and Hall/CRC.

See Also

rmh, rStrauss, rStraussHard, rDiggleGratton. rDGS, rPenttinen.

For fitting the model, see Hardcore.

Examples

X <- rHardcore(0.05,1.5,square(50))

---

Random Perturbation of Line Segment Pattern

Description

Randomly pertubs a spatial pattern of line segments by applying independent random displacements to the segment endpoints.

Usage

## S3 method for class 'psp'
rjitter(X, radius, ..., clip=TRUE, nsim=1, drop=TRUE)

Arguments

X	A point pattern on a linear network (object of class "psp").
radius	Scale of perturbations. A positive numerical value. Each point will be displaced by a random distance, with maximum displacement equal to this value.
...	Ignored.
clip	Logical value specifying what to do if segments cross the boundary of the window. See Details.
nsim	Number of simulated realisations to be generated.
drop	Logical. If nsim=1 and drop=TRUE (the default), the result will be a spatial pattern of line segments (class "psp") rather than a list of length 1 containing this pattern.

Details

The function rjitter is generic. This function is the method for the class "psp" of line segment patterns.

Each of the endpoints of each segment in X will be subjected to an independent random displacement. The displacement vectors are uniformly distributed in a circle of radius radius.

If clip=TRUE (the default), segment endpoints are permitted to move to locations slightly outside the window of X, and the resulting segments will be clipped to the window. If clip=FALSE, segment endpoints are conditioned to fall inside the window.

If nsim=1 and drop=TRUE, the result is another spatial pattern of line segments (object of class "psp"). Otherwise, the result is a list of nsim line segment patterns.

Value

A spatial pattern of line segments (object of class "psp") or a list of such patterns.

Author(s)

Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

See Also

rjitter for point patterns in two dimensions.

Examples

E <- edges(letterR)
   Window(E) <- owin(c(1.9, 4.1), c(0.5, 3.5))
   plot(rjitter(E, 0.1))

---

Theoretical Distribution of Nearest Neighbour Distance

Description

Density, distribution function, quantile function and random generation for the random distance to the $k$th nearest neighbour in a Poisson point process in $d$ dimensions.

Usage

dknn(x, k = 1, d = 2, lambda = 1)
pknn(q, k = 1, d = 2, lambda = 1)
qknn(p, k = 1, d = 2, lambda = 1)
rknn(n, k = 1, d = 2, lambda = 1)

Arguments

x, q	vector of quantiles.
p	vector of probabilities.
n	number of observations to be generated.
k	order of neighbour.
d	dimension of space.
lambda	intensity of Poisson point process.

Details

In a Poisson point process in $d$-dimensional space, let the random variable $R$ be the distance from a fixed point to the $k$-th nearest random point, or the distance from a random point to the $k$-th nearest other random point.

Then $R^d$ has a Gamma distribution with shape parameter $k$ and rate $\lambda * \alpha$ where $\alpha$ is a constant (equal to the volume of the unit ball in $d$-dimensional space). See e.g. Cressie (1991, page 61).

These functions support calculation and simulation for the distribution of $R$.

Value

A numeric vector: dknn returns the probability density, pknn returns cumulative probabilities (distribution function), qknn returns quantiles, and rknn generates random deviates.

Author(s)

Adrian Baddeley [email protected]

and Rolf Turner [email protected]

References

Cressie, N.A.C. (1991) Statistics for spatial data. John Wiley and Sons, 1991.

Examples

x <- seq(0, 5, length=20)
  densities <- dknn(x, k=3, d=2)
  cdfvalues <- pknn(x, k=3, d=2)
  randomvalues <- rknn(100, k=3, d=2)
  deciles <- qknn((1:9)/10, k=3, d=2)

---

Random Re-Labelling of Point Pattern

Description

Randomly allocates marks to a point pattern, or permutes the existing marks, or resamples from the existing marks.

Usage

rlabel(X, labels=marks(X), permute=TRUE, group=NULL, ..., nsim=1, drop=TRUE)

Arguments

X	Point pattern (object of class "ppp", "lpp", "pp3" or "ppx") or line segment pattern (object of class "psp").
labels	Vector of values from which the new marks will be drawn at random. Defaults to the vector of existing marks.
permute	Logical value indicating whether to generate new marks by randomly permuting labels or by drawing a random sample with replacement.
group	Optional. A factor, or other data dividing the points into groups. Random relabelling will be performed separately within each group. See Details.
...	Additional arguments passed to cut.ppp to determine the grouping factor, when group is given.
nsim	Number of simulated realisations to be generated.
drop	Logical. If nsim=1 and drop=TRUE (the default), the result will be a point pattern, rather than a list containing a point pattern.

Details

This very simple function allocates random marks to an existing point pattern X. It is useful for hypothesis testing purposes. (The function can also be applied to line segment patterns.)

In the simplest case, the command rlabel(X) yields a point pattern obtained from X by randomly permuting the marks of the points.

If permute=TRUE, then labels should be a vector of length equal to the number of points in X. The result of rlabel will be a point pattern with locations given by X and marks given by a random permutation of labels (i.e. a random sample without replacement).

If permute=FALSE, then labels may be a vector of any length. The result of rlabel will be a point pattern with locations given by X and marks given by a random sample from labels (with replacement).

The argument group specifies that the points are divided into several different groups, and that the random labelling shall be performed separately on each group. The arguments group and ... are passed to cut.ppp to determine the grouping. Thus group could be a factor, or the name of a column of marks in X, or a tessellation, or a factor-valued pixel image, etc.

Value

If nsim = 1 and drop=TRUE, a marked point pattern (of the same class as X). If nsim > 1, a list of point patterns.

Author(s)

Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

See Also

marks<- to assign arbitrary marks.

Examples

amacrine

   # Randomly permute the marks "on" and "off"
   # Result always has 142 "off" and 152 "on"
   Y <- rlabel(amacrine)

   # randomly allocate marks "on" and "off"
   # with probabilities p(off) = 0.48, p(on) = 0.52
   Y <- rlabel(amacrine, permute=FALSE)

   # randomly allocate marks "A" and "B" with equal probability
   Y <- rlabel(cells, labels=factor(c("A", "B")), permute=FALSE)

   # divide the window into tiles and 
   # randomly permute the marks within each tile
   Z <- rlabel(amacrine, group=quadrats(Window(amacrine), 4, 3))

---

Simulate Log-Gaussian Cox Process

Description

Generate a random point pattern, a realisation of the log-Gaussian Cox process.

Usage

rLGCP(model=c("exponential", "gauss", "stable", "gencauchy", "matern"),
       mu = 0, param = NULL,
       ...,
       win=NULL, saveLambda=TRUE, nsim=1, drop=TRUE)

Arguments

model	character string (partially matched) giving the name of a covariance model for the Gaussian random field.
mu	mean function of the Gaussian random field. Either a single number, a function(x,y, ...) or a pixel image (object of class "im").
param	List of parameters for the covariance. Standard arguments are var and scale.
...	Additional parameters for the covariance, or arguments passed to as.mask to determine the pixel resolution.
win	Window in which to simulate the pattern. An object of class "owin".
saveLambda	Logical. If TRUE (the default) then the simulated random intensity will also be saved, and returns as an attribute of the point pattern.
nsim	Number of simulated realisations to be generated.
drop	Logical. If nsim=1 and drop=TRUE (the default), the result will be a point pattern, rather than a list containing a point pattern.

Details

This function generates a realisation of a log-Gaussian Cox process (LGCP). This is a Cox point process in which the logarithm of the random intensity is a Gaussian random field with mean function $\mu$ and covariance function $c(r)$. Conditional on the random intensity, the point process is a Poisson process with this intensity.

The string model specifies the covariance function of the Gaussian random field, and the parameters of the covariance are determined by param and ....

All models recognise the parameters var for the variance at distance zero, and scale for the scale parameter. Some models require additional parameters which are listed below.

The available models are as follows:

model="exponential":

the exponential covariance function

$C(r) = \sigma^2 \exp(-r/h)$

where $\sigma^2$ is the variance parameter var, and $h$ is the scale parameter scale.

model="gauss":

the Gaussian covariance function

$C(r) = \sigma^2 \exp(-(r/h)^2)$

where $\sigma^2$ is the variance parameter var, and $h$ is the scale parameter scale.

model="stable":

the stable covariance function

$C(r) = \sigma^2 \exp(-(r/h)^\alpha)$

where $\sigma^2$ is the variance parameter var, $h$ is the scale parameter scale, and $\alpha$ is the shape parameter alpha. The parameter alpha must be given, either as a stand-alone argument, or as an entry in the list param.

model="gencauchy":

the generalised Cauchy covariance function

$C(r) = \sigma^2 (1 + (x/h)^\alpha)^{-\beta/\alpha}$

where $\sigma^2$ is the variance parameter var, $h$ is the scale parameter scale, and $\alpha$ and $\beta$ are the shape parameters alpha and beta. The parameters alpha and beta must be given, either as stand-alone arguments, or as entries in the list param.

model="matern":

the Whittle-Matern covariance function

$C(r) = \sigma^2 \frac{1}{2^{\nu-1} \Gamma(\nu)} (\sqrt{2 \nu} \, r/h)^\nu K_\nu(\sqrt{2\nu}\, r/h)$

where $\sigma^2$ is the variance parameter var, $h$ is the scale parameter scale, and $\nu$ is the shape parameter nu. The parameter nu must be given, either as a stand-alone argument, or as an entry in the list param.

The algorithm uses the circulant embedding technique to generate values of a Gaussian random field, with the specified mean function mu and the covariance specified by the arguments model and param, on the points of a regular grid. The exponential of this random field is taken as the intensity of a Poisson point process, and a realisation of the Poisson process is then generated by the function rpoispp in the spatstat.random package.

If the simulation window win is missing or NULL, then it defaults to Window(mu) if mu is a pixel image, and it defaults to the unit square otherwise.

The LGCP model can be fitted to data using kppm.

Value

A point pattern (object of class "ppp") or a list of point patterns.

Additionally, the simulated intensity function for each point pattern is returned as an attribute "Lambda" of the point pattern, if saveLambda=TRUE.

Warning: new implementation

The simulation algorithm for rLGCP has been completely re-written in spatstat.random version 3.2-0 to avoid depending on the package RandomFields which is now defunct (and is sadly missed).

It is no longer possible to replicate results that were obtained using rLGCP in previous versions of spatstat.random.

The current code is a new implementation and should be considered vulnerable to new bugs.

Author(s)

Abdollah Jalilian and Rasmus Waagepetersen. Modified by Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

References

Moller, J., Syversveen, A. and Waagepetersen, R. (1998) Log Gaussian Cox Processes. Scandinavian Journal of Statistics 25, 451–482.

See Also

rpoispp, rMatClust, rGaussPoisson, rNeymanScott.

For fitting the model, see kppm, lgcp.estK.

Examples

online <- interactive()

    # homogeneous LGCP with exponential covariance function
    X <- rLGCP("exp", 3, var=0.2, scale=.1)

    # inhomogeneous LGCP with Gaussian covariance function
    m <- as.im(function(x, y){5 - 1.5 * (x - 0.5)^2 + 2 * (y - 0.5)^2}, W=owin())
    X <- rLGCP("gauss", m, var=0.15, scale =0.1)

    if(online) {
     plot(attr(X, "Lambda"))
     points(X)
    }

    # inhomogeneous LGCP with Matern covariance function
    X <- rLGCP("matern", function(x, y){ 1 - 0.4 * x},
               var=2, scale=0.7, nu=0.5,
               win = owin(c(0, 10), c(0, 10)))
    if(online) plot(X)

---

Simulate Matern Cluster Process

Description

Generate a random point pattern, a simulated realisation of the Matern Cluster Process.

Usage

rMatClust(kappa, scale, mu, win = square(1),
           nsim=1, drop=TRUE, ...,
           n.cond=NULL, w.cond=NULL,
           algorithm=c("BKBC", "naive"),
           nonempty=TRUE, 
           poisthresh=1e-6, saveparents=FALSE, saveLambda=FALSE,
           kappamax=NULL, mumax=NULL)

Arguments

kappa	Intensity of the Poisson process of cluster centres. A single positive number, a function, or a pixel image.
scale	Radius of the clusters. A single positive number.
mu	Mean number of points per cluster (a single positive number) or reference intensity for the cluster points (a function or a pixel image).
win	Window in which to simulate the pattern. An object of class "owin" or something acceptable to as.owin.
nsim	Number of simulated realisations to be generated.
drop	Logical. If nsim=1 and drop=TRUE (the default), the result will be a point pattern, rather than a list containing a point pattern.
...	Passed to clusterfield to control the image resolution when saveLambda=TRUE.
n.cond	Optional. Integer specifying a fixed number of points. See the section on Conditional Simulation.
w.cond	Optional. Conditioning region. A window (object of class "owin") specifying the region which must contain exactly n.cond points. See the section on Conditional Simulation.
algorithm	String (partially matched) specifying the simulation algorithm. See Details.
nonempty	Logical. If TRUE (the default), a more efficient algorithm is used, in which parents are generated conditionally on having at least one offspring point in the window. If FALSE, parents are generated even if they have no offspring in the window. The default is recommended unless you need to simulate all the parent points for some other purpose.
poisthresh	Numerical threshold below which the model will be treated as a Poisson process. See Details.
saveparents	Logical value indicating whether to save the locations of the parent points as an attribute.
saveLambda	Logical. If TRUE then the random intensity corresponding to the simulated parent points will also be calculated and saved, and returns as an attribute of the point pattern.
kappamax	Optional. Numerical value which is an upper bound for the values of kappa, when kappa is a pixel image or a function.
mumax	Optional. Numerical value which is an upper bound for the values of mu, when mu is a pixel image or a function.

Details

This algorithm generates a realisation of Matern's cluster process, a special case of the Neyman-Scott process, inside the window win.

In the simplest case, where kappa and mu are single numbers, the cluster process is formed by first generating a uniform Poisson point process of “parent” points with intensity kappa. Then each parent point is replaced by a random cluster of “offspring” points, the number of points per cluster being Poisson (mu) distributed, and their positions being placed and uniformly inside a disc of radius scale centred on the parent point. The resulting point pattern is a realisation of the classical “stationary Matern cluster process”. This point process has intensity kappa * mu.

The algorithm can also generate spatially inhomogeneous versions of the Matern cluster process:

• The parent points can be spatially inhomogeneous. If the argument kappa is a function(x,y) or a pixel image (object of class "im"), then it is taken as specifying the intensity function of an inhomogeneous Poisson process that generates the parent points.

• The offspring points can be inhomogeneous. If the argument mu is a function(x,y) or a pixel image (object of class "im"), then it is interpreted as the reference density for offspring points, in the sense of Waagepetersen (2007). For a given parent point, the offspring constitute a Poisson process with intensity function equal to mu/(pi * scale^2) inside the disc of radius scale centred on the parent point, and zero intensity outside this disc. Equivalently we first generate, for each parent point, a Poisson ($M$) random number of offspring (where $M$ is the maximum value of mu) placed independently and uniformly in the disc of radius scale centred on the parent location, and then randomly thin the offspring points, with retention probability mu/M.

• Both the parent points and the offspring points can be inhomogeneous, as described above.

The intensity of the Matern cluster process is kappa * mu if either kappa or mu is a single number. In the general case the intensity is an integral involving kappa, mu and scale.

Value

A point pattern (an object of class "ppp") if nsim=1, or a list of point patterns if nsim > 1.

Additionally, some intermediate results of the simulation are returned as attributes of this point pattern (see rNeymanScott). Furthermore, the simulated intensity function is returned as an attribute "Lambda", if saveLambda=TRUE.

Simulation Algorithm

Two simulation algorithms are implemented.

• The naive algorithm generates the cluster process by directly following the description given above. First the window win is expanded by a distance equal to scale. Then the parent points are generated in the expanded window according to a Poisson process with intensity kappa. Then each parent point is replaced by a finite cluster of offspring points as described above. The naive algorithm is used if algorithm="naive" or if nonempty=FALSE.

• The BKBC algorithm, proposed by Baddeley and Chang (2023), is a modification of the algorithm of Brix and Kendall (2002). Parents are generated in the infinite plane, subject to the condition that they have at least one offspring point inside the window win. The BKBC algorithm is used when algorithm="BKBC" (the default) and nonempty=TRUE (the default).

The naive algorithm becomes very slow when scale is large, while the BKBC algorithm is uniformly fast (Baddeley and Chang, 2023).

If saveparents=TRUE, then the simulated point pattern will have an attribute "parents" containing the coordinates of the parent points, and an attribute "parentid" mapping each offspring point to its parent.

If nonempty=TRUE (the default), then parents are generated subject to the condition that they have at least one offspring point in the window win. nonempty=FALSE, then parents without offspring will be included; this option is not available in the BKBC algorithm.

Note that if kappa is a pixel image, its domain must be larger than the window win. This is because an offspring point inside win could have its parent point lying outside win. In order to allow this, the naive simulation algorithm first expands the original window win by a distance equal to scale and generates the Poisson process of parent points on this larger window. If kappa is a pixel image, its domain must contain this larger window.

If the pair correlation function of the model is very close to that of a Poisson process, with maximum deviation less than poisthresh, then the model is approximately a Poisson process. This is detected by the naive algorithm which then simulates a Poisson process with intensity kappa * mu, using rpoispp. This avoids computations that would otherwise require huge amounts of memory.

Fitting cluster models to data

The Matern cluster process model with homogeneous parents (i.e. where kappa is a single number) where the offspring are either homogeneous or inhomogeneous (mu is a single number, a function or pixel image) can be fitted to point pattern data using kppm, or fitted to the inhomogeneous $K$ function using matclust.estK or matclust.estpcf.

Currently spatstat does not support fitting the Matern cluster process model with inhomogeneous parents.

A fitted Matern cluster process model can be simulated automatically using simulate.kppm (which invokes rMatClust to perform the simulation).

Conditional Simulation

If n.cond is specified, it should be a single integer. Simulation will be conditional on the event that the pattern contains exactly n.cond points (or contains exactly n.cond points inside the region w.cond if it is given).

Conditional simulation uses the rejection algorithm described in Section 6.2 of Moller, Syversveen and Waagepetersen (1998). There is a maximum number of proposals which will be attempted. Consequently the return value may contain fewer than nsim point patterns.

Author(s)

Adrian Baddeley [email protected], Ya-Mei Chang [email protected] and Rolf Turner [email protected].

References

Baddeley, A. and Chang, Y.-M. (2023) Robust algorithms for simulating cluster point processes. Journal of Statistical Computation and Simulation. In Press. DOI 10.1080/00949655.2023.2166045.

Brix, A. and Kendall, W.S. (2002) Simulation of cluster point processes without edge effects. Advances in Applied Probability 34, 267–280.

Matern, B. (1960) Spatial Variation. Meddelanden fraan Statens Skogsforskningsinstitut, volume 59, number 5. Statens Skogsforskningsinstitut, Sweden.

Matern, B. (1986) Spatial Variation. Lecture Notes in Statistics 36, Springer-Verlag, New York.

Moller, J., Syversveen, A. and Waagepetersen, R. (1998) Log Gaussian Cox Processes. Scandinavian Journal of Statistics 25, 451–482.

Waagepetersen, R. (2007) An estimating function approach to inference for inhomogeneous Neyman-Scott processes. Biometrics 63, 252–258.

See Also

rpoispp, rThomas, rCauchy, rVarGamma, rNeymanScott, rGaussPoisson.

For fitting the model, see kppm, clusterfit.

Examples

# homogeneous
 X <- rMatClust(10, 0.05, 4)
 # inhomogeneous
 ff <- function(x,y){ 4 * exp(2 * abs(x) - 1) }
 Z <- as.im(ff, owin())
 Y <- rMatClust(10, 0.05, Z)
 YY <- rMatClust(ff, 0.05, 3)

---

Simulate Matern Model I

Description

Generate a random point pattern, a simulated realisation of the Matern Model I inhibition process model.

Usage

rMaternI(kappa, r, win = owin(c(0,1),c(0,1)), stationary=TRUE, ...,
          nsim=1, drop=TRUE)

Arguments

kappa	Intensity of the Poisson process of proposal points. A single positive number.
r	Inhibition distance.
win	Window in which to simulate the pattern. An object of class "owin" or something acceptable to as.owin. Alternatively a higher-dimensional box of class "box3" or "boxx".
stationary	Logical. Whether to start with a stationary process of proposal points (stationary=TRUE) or to generate the proposal points only inside the window (stationary=FALSE).
...	Ignored.
nsim	Number of simulated realisations to be generated.
drop	Logical. If nsim=1 and drop=TRUE (the default), the result will be a point pattern, rather than a list containing a point pattern.

Details

This algorithm generates one or more realisations of Matern's Model I inhibition process inside the window win.

The process is constructed by first generating a uniform Poisson point process of “proposal” points with intensity kappa. If stationary = TRUE (the default), the proposal points are generated in a window larger than win that effectively means the proposals are stationary. If stationary=FALSE then the proposal points are only generated inside the window win.

A proposal point is then deleted if it lies within r units' distance of another proposal point. Otherwise it is retained.

The retained points constitute Matern's Model I.

Value

A point pattern if nsim=1, or a list of point patterns if nsim > 1. Each point pattern is normally an object of class "ppp", but may be of class "pp3" or "ppx" depending on the window.

Author(s)

Ute Hahn, Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

See Also

rpoispp, rMatClust

Examples

X <- rMaternI(20, 0.05)
 Y <- rMaternI(20, 0.05, stationary=FALSE)

---

Simulate Matern Model II

Description

Generate a random point pattern, a simulated realisation of the Matern Model II inhibition process.

Usage

rMaternII(kappa, r, win = owin(c(0,1),c(0,1)), stationary=TRUE, ...,
           nsim=1, drop=TRUE)

Arguments

kappa	Intensity of the Poisson process of proposal points. A single positive number.
r	Inhibition distance.
win	Window in which to simulate the pattern. An object of class "owin" or something acceptable to as.owin. Alternatively a higher-dimensional box of class "box3" or "boxx".
stationary	Logical. Whether to start with a stationary process of proposal points (stationary=TRUE) or to generate the proposal points only inside the window (stationary=FALSE).
...	Ignored.
nsim	Number of simulated realisations to be generated.
drop	Logical. If nsim=1 and drop=TRUE (the default), the result will be a point pattern, rather than a list containing a point pattern.

Details

This algorithm generates one or more realisations of Matern's Model II inhibition process inside the window win.

The process is constructed by first generating a uniform Poisson point process of “proposal” points with intensity kappa. If stationary = TRUE (the default), the proposal points are generated in a window larger than win that effectively means the proposals are stationary. If stationary=FALSE then the proposal points are only generated inside the window win.

Then each proposal point is marked by an “arrival time”, a number uniformly distributed in $[0,1]$ independently of other variables.

A proposal point is deleted if it lies within r units' distance of another proposal point that has an earlier arrival time. Otherwise it is retained. The retained points constitute Matern's Model II.

The difference between Matern's Model I and II is the italicised statement above. Model II has a higher intensity for the same parameter values.

Value

A point pattern if nsim=1, or a list of point patterns if nsim > 1. Each point pattern is normally an object of class "ppp", but may be of class "pp3" or "ppx" depending on the window.

Author(s)

Ute Hahn, Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

See Also

rpoispp, rMatClust, rMaternI

Examples

X <- rMaternII(20, 0.05)
 Y <- rMaternII(20, 0.05, stationary=FALSE)

---

Simulate point patterns using the Metropolis-Hastings algorithm.

Description

Generic function for running the Metropolis-Hastings algorithm to produce simulated realisations of a point process model.

Usage

rmh(model, ...)

Arguments

model	The point process model to be simulated.
...	Further arguments controlling the simulation.

Details

The Metropolis-Hastings algorithm can be used to generate simulated realisations from a wide range of spatial point processes. For caveats, see below.

The function rmh is generic; it has methods rmh.ppm (for objects of class "ppm") and rmh.default (the default). The actual implementation of the Metropolis-Hastings algorithm is contained in rmh.default. For details of its use, see rmh.ppm or rmh.default.

[If the model is a Poisson process, then Metropolis-Hastings is not used; the Poisson model is generated directly using rpoispp or rmpoispp.]

In brief, the Metropolis-Hastings algorithm is a Markov Chain, whose states are spatial point patterns, and whose limiting distribution is the desired point process. After running the algorithm for a very large number of iterations, we may regard the state of the algorithm as a realisation from the desired point process.

However, there are difficulties in deciding whether the algorithm has run for “long enough”. The convergence of the algorithm may indeed be extremely slow. No guarantees of convergence are given!

While it is fashionable to decry the Metropolis-Hastings algorithm for its poor convergence and other properties, it has the advantage of being easy to implement for a wide range of models.

Value

A point pattern, in the form of an object of class "ppp". See rmh.default for details.

Warning

As of version 1.22-1 of spatstat a subtle change was made to rmh.default(). We had noticed that the results produced were sometimes not “scalable” in that two models, differing in effect only by the units in which distances are measured and starting from the same seed, gave different results. This was traced to an idiosyncracy of floating point arithmetic. The code of rmh.default() has been changed so that the results produced by rmh are now scalable. The downside of this is that code which users previously ran may now give results which are different from what they formerly were.

In order to recover former behaviour (so that previous results can be reproduced) set spatstat.options(scalable=FALSE). See the last example in the help for rmh.default.

Author(s)

Adrian Baddeley [email protected]

and Rolf Turner [email protected]

See Also

rmh.default

Examples

# See examples in rmh.default and rmh.ppm

---

Simulate Point Process Models using the Metropolis-Hastings Algorithm.

Description

Generates a random point pattern, simulated from a chosen point process model, using the Metropolis-Hastings algorithm.

Usage

## Default S3 method:
rmh(model, start=NULL,
   control=default.rmhcontrol(model),
   ...,
   nsim=1, drop=TRUE, saveinfo=TRUE,
   verbose=TRUE, snoop=FALSE)

Arguments

model	Data specifying the point process model that is to be simulated.
start	Data determining the initial state of the algorithm.
control	Data controlling the iterative behaviour and termination of the algorithm.
...	Further arguments passed to rmhcontrol or to trend functions in model.
nsim	Number of simulated point patterns that should be generated.
drop	Logical. If nsim=1 and drop=TRUE (the default), the result will be a point pattern, rather than a list containing a single point pattern.
saveinfo	Logical value indicating whether to save auxiliary information.
verbose	Logical value indicating whether to print progress reports.
snoop	Logical. If TRUE, activate the visual debugger.

Details

This function generates simulated realisations from any of a range of spatial point processes, using the Metropolis-Hastings algorithm. It is the default method for the generic function rmh.

This function executes a Metropolis-Hastings algorithm with birth, death and shift proposals as described in Geyer and Moller (1994).

The argument model specifies the point process model to be simulated. It is either a list, or an object of class "rmhmodel", with the following components:

cif

A character string specifying the choice of interpoint interaction for the point process.

par

Parameter values for the conditional intensity function.

w

(Optional) window in which the pattern is to be generated. An object of class "owin", or data acceptable to as.owin.

trend

Data specifying the spatial trend in the model, if it has a trend. This may be a function, a pixel image (of class "im"), (or a list of functions or images if the model is multitype).

If the trend is a function or functions, any auxiliary arguments ... to rmh.default will be passed to these functions, which should be of the form function(x, y, ...).

types

List of possible types, for a multitype point process.

For full details of these parameters, see rmhmodel.default.

The argument start determines the initial state of the Metropolis-Hastings algorithm. It is either NULL, or an object of class "rmhstart", or a list with the following components:

n.start

Number of points in the initial point pattern. A single integer, or a vector of integers giving the numbers of points of each type in a multitype point pattern. Incompatible with x.start.

x.start

Initial point pattern configuration. Incompatible with n.start.

x.start may be a point pattern (an object of class "ppp"), or data which can be coerced to this class by as.ppp, or an object with components x and y, or a two-column matrix. In the last two cases, the window for the pattern is determined by model$w. In the first two cases, if model$w is also present, then the final simulated pattern will be clipped to the window model$w.

For full details of these parameters, see rmhstart.

The third argument control controls the simulation procedure (including conditional simulation), iterative behaviour, and termination of the Metropolis-Hastings algorithm. It is either NULL, or a list, or an object of class "rmhcontrol", with components:

p

The probability of proposing a “shift” (as opposed to a birth or death) in the Metropolis-Hastings algorithm.

q

The conditional probability of proposing a death (rather than a birth) given that birth/death has been chosen over shift.

nrep

The number of repetitions or iterations to be made by the Metropolis-Hastings algorithm. It should be large.

expand

Either a numerical expansion factor, or a window (object of class "owin"). Indicates that the process is to be simulated on a larger domain than the original data window w, then clipped to w when the algorithm has finished.

The default is to expand the simulation window if the model is stationary and non-Poisson (i.e. it has no trend and the interaction is not Poisson) and not to expand in all other cases.

If the model has a trend, then in order for expansion to be feasible, the trend must be given either as a function, or an image whose bounding box is large enough to contain the expanded window.

periodic

A logical scalar; if periodic is TRUE we simulate a process on the torus formed by identifying opposite edges of a rectangular window.

ptypes

A vector of probabilities (summing to 1) to be used in assigning a random type to a new point.

fixall

A logical scalar specifying whether to condition on the number of points of each type.

nverb

An integer specifying how often “progress reports” (which consist simply of the number of repetitions completed) should be printed out. If nverb is left at 0, the default, the simulation proceeds silently.

x.cond

If this argument is present, then conditional simulation will be performed, and x.cond specifies the conditioning points and the type of conditioning.

nsave,nburn

If these values are specified, then intermediate states of the simulation algorithm will be saved every nsave iterations, after an initial burn-in period of nburn iterations.

track

Logical flag indicating whether to save the transition history of the simulations.

For full details of these parameters, see rmhcontrol. The control parameters can also be given in the ... arguments.

Value

A point pattern (an object of class "ppp", see ppp.object) or a list of point patterns.

The returned value has an attribute info containing modified versions of the arguments model, start, and control which together specify the exact simulation procedure. The info attribute can be printed (and is printed automatically by summary.ppp). For computational efficiency, the info attribute can be omitted by setting saveinfo=FALSE.

The value of .Random.seed at the start of the simulations is also saved and returned as an attribute seed.

If the argument track=TRUE was given (see rmhcontrol), the transition history of the algorithm is saved, and returned as an attribute history. The transition history is a data frame containing a factor proposaltype identifying the proposal type (Birth, Death or Shift) and a logical vector accepted indicating whether the proposal was accepted. The data frame also has columns numerator, denominator which give the numerator and denominator of the Hastings ratio for the proposal.

If the argument nsave was given (see rmhcontrol), the return value has an attribute saved which is a list of point patterns, containing the intermediate states of the algorithm.

Conditional Simulation

There are several kinds of conditional simulation.

• Simulation conditional upon the number of points, that is, holding the number of points fixed. To do this, set control$p (the probability of a shift) equal to 1. The number of points is then determined by the starting state, which may be specified either by setting start$n.start to be a scalar, or by setting the initial pattern start$x.start.

• In the case of multitype processes, it is possible to simulate the model conditionally upon the number of points of each type, i.e. holding the number of points of each type to be fixed. To do this, set control$p equal to 1 and control$fixall to be TRUE. The number of points is then determined by the starting state, which may be specified either by setting start$n.start to be an integer vector, or by setting the initial pattern start$x.start.

• Simulation conditional on the configuration observed in a sub-window, that is, requiring that, inside a specified sub-window $V$, the simulated pattern should agree with a specified point pattern $y$.To do this, set control$x.cond to equal the specified point pattern $y$, making sure that it is an object of class "ppp" and that the window Window(control$x.cond) is the conditioning window $V$.

• Simulation conditional on the presence of specified points, that is, requiring that the simulated pattern should include a specified set of points. This is simulation from the Palm distribution of the point process given a pattern $y$. To do this, set control$x.cond to be a data.frame containing the coordinates (and marks, if appropriate) of the specified points.

For further information, see rmhcontrol.

Note that, when we simulate conditionally on the number of points, or conditionally on the number of points of each type, no expansion of the window is possible.

Visual Debugger

If snoop = TRUE, an interactive debugger is activated. On the current plot device, the debugger displays the current state of the Metropolis-Hastings algorithm together with the proposed transition to the next state. Clicking on this graphical display (using the left mouse button) will re-centre the display at the clicked location. Surrounding this graphical display is an array of boxes representing different actions. Clicking on one of the action boxes (using the left mouse button) will cause the action to be performed. Debugger actions include:

• Zooming in or out

• Panning (shifting the field of view) left, right, up or down

• Jumping to the next iteration

• Skipping 10, 100, 1000, 10000 or 100000 iterations

• Jumping to the next Birth proposal (etc)

• Changing the fate of the proposal (i.e. changing whether the proposal is accepted or rejected)

• Dumping the current state and proposal to a file

• Printing detailed information at the terminal

• Exiting the debugger (so that the simulation algorithm continues without further interruption).

Right-clicking the mouse will also cause the debugger to exit.

Warnings

There is never a guarantee that the Metropolis-Hastings algorithm has converged to its limiting distribution.

If start$x.start is specified then expand is set equal to 1 and simulation takes place in Window(x.start). Any specified value for expand is simply ignored.

The presence of both a component w of model and a non-null value for Window(x.start) makes sense ONLY if w is contained in Window(x.start).

For multitype processes make sure that, even if there is to be no trend corresponding to a particular type, there is still a component (a NULL component) for that type, in the list.

Other models

In theory, any finite point process model can be simulated using the Metropolis-Hastings algorithm, provided the conditional intensity is uniformly bounded.

In practice, the list of point process models that can be simulated using rmh.default is limited to those that have been implemented in the package's internal C code. More options will be added in the future.

Note that the lookup conditional intensity function permits the simulation (in theory, to any desired degree of approximation) of any pairwise interaction process for which the interaction depends only on the distance between the pair of points.

Reproducible simulations

If the user wants the simulation to be exactly reproducible (e.g. for a figure in a journal article, where it is useful to have the figure consistent from draft to draft) then the state of the random number generator should be set before calling rmh.default. This can be done either by calling set.seed or by assigning a value to .Random.seed. In the examples below, we use set.seed.

If a simulation has been performed and the user now wants to repeat it exactly, the random seed should be extracted from the simulated point pattern X by seed <- attr(x, "seed"), then assigned to the system random nunber state by .Random.seed <- seed before calling rmh.default.

Author(s)

Adrian Baddeley [email protected] and Rolf Turner [email protected]

References

Baddeley, A. and Turner, R. (2000) Practical maximum pseudolikelihood for spatial point patterns. Australian and New Zealand Journal of Statistics 42, 283 – 322.

Diggle, P. J. (2003) Statistical Analysis of Spatial Point Patterns (2nd ed.) Arnold, London.

Diggle, P.J. and Gratton, R.J. (1984) Monte Carlo methods of inference for implicit statistical models. Journal of the Royal Statistical Society, series B 46, 193 – 212.

Diggle, P.J., Gates, D.J., and Stibbard, A. (1987) A nonparametric estimator for pairwise-interaction point processes. Biometrika 74, 763 – 770.

Geyer, C.J. and Moller, J. (1994) Simulation procedures and likelihood inference for spatial point processes. Scandinavian Journal of Statistics 21, 359–373.

Geyer, C.J. (1999) Likelihood Inference for Spatial Point Processes. Chapter 3 in O.E. Barndorff-Nielsen, W.S. Kendall and M.N.M. Van Lieshout (eds) Stochastic Geometry: Likelihood and Computation, Chapman and Hall / CRC, Monographs on Statistics and Applied Probability, number 80. Pages 79–140.

See Also

rmh, rmh.ppm, rStrauss, ppp, ppm

Interactions: AreaInter, BadGey, DiggleGatesStibbard, DiggleGratton, Fiksel, Geyer, Hardcore, Hybrid, LennardJones, MultiStrauss, MultiStraussHard, PairPiece, Penttinen, Poisson, Softcore, Strauss, StraussHard and Triplets.

Examples

if(interactive()) {
     nr   <- 1e5
     nv  <- 5000
     ns <- 200
   } else {
     nr  <- 20
     nv <- 5
     ns <- 20
     oldopt <- spatstat.options()
     spatstat.options(expand=1.05)
   }
   set.seed(961018)
   
   # Strauss process.
   mod01 <- list(cif="strauss",par=list(beta=2,gamma=0.2,r=0.7),
                 w=c(0,10,0,10))
   X1.strauss <- rmh(model=mod01,start=list(n.start=ns),
                     control=list(nrep=nr,nverb=nv))

   if(interactive()) plot(X1.strauss)
   
   # Strauss process, conditioning on n = 42:
   X2.strauss <- rmh(model=mod01,start=list(n.start=42),
                     control=list(p=1,nrep=nr,nverb=nv))

   # Tracking algorithm progress:
   # (a) saving intermediate states:
   X <- rmh(model=mod01,start=list(n.start=ns),
            control=list(nrep=nr, nsave=nr/5, nburn=nr/2))
   Saved <- attr(X, "saved")
   plot(Saved)

   # (b) inspecting transition history:
   X <- rmh(model=mod01,start=list(n.start=ns),
            control=list(nrep=nr, track=TRUE))
   History <- attr(X, "history")
   head(History)

   # Hard core process:
   mod02 <- list(cif="hardcore",par=list(beta=2,hc=0.7),w=c(0,10,0,10))
   X3.hardcore <- rmh(model=mod02,start=list(n.start=ns),
                     control=list(nrep=nr,nverb=nv))
   
   if(interactive()) plot(X3.hardcore)

   # Strauss process equal to pure hardcore:
   mod02s <- list(cif="strauss",par=list(beta=2,gamma=0,r=0.7),w=c(0,10,0,10))
   X3.strauss <- rmh(model=mod02s,start=list(n.start=ns),
                     control=list(nrep=nr,nverb=nv))
   
   # Strauss process in a polygonal window.
   x     <- c(0.55,0.68,0.75,0.58,0.39,0.37,0.19,0.26,0.42)
   y     <- c(0.20,0.27,0.68,0.99,0.80,0.61,0.45,0.28,0.33)
   mod03 <- list(cif="strauss",par=list(beta=2000,gamma=0.6,r=0.07),
                w=owin(poly=list(x=x,y=y)))
   X4.strauss <- rmh(model=mod03,start=list(n.start=ns),
                     control=list(nrep=nr,nverb=nv))
   if(interactive()) plot(X4.strauss)
   
   # Strauss process in a polygonal window, conditioning on n = 80.
   X5.strauss <- rmh(model=mod03,start=list(n.start=ns),
                     control=list(p=1,nrep=nr,nverb=nv))
   
   # Strauss process, starting off from X4.strauss, but with the
   # polygonal window replace by a rectangular one.  At the end,
   # the generated pattern is clipped to the original polygonal window.
   xxx <- X4.strauss
   Window(xxx) <- as.owin(c(0,1,0,1))
   X6.strauss <- rmh(model=mod03,start=list(x.start=xxx),
                     control=list(nrep=nr,nverb=nv))
   
   # Strauss with hardcore:
   mod04 <- list(cif="straush",par=list(beta=2,gamma=0.2,r=0.7,hc=0.3),
                w=c(0,10,0,10))
   X1.straush <- rmh(model=mod04,start=list(n.start=ns),
                     control=list(nrep=nr,nverb=nv))
   
   # Another Strauss with hardcore (with a perhaps surprising result):
   mod05 <- list(cif="straush",par=list(beta=80,gamma=0.36,r=45,hc=2.5),
                w=c(0,250,0,250))
   X2.straush <- rmh(model=mod05,start=list(n.start=ns),
                     control=list(nrep=nr,nverb=nv))
   
   # Pure hardcore (identical to X3.strauss).
   mod06 <- list(cif="straush",par=list(beta=2,gamma=1,r=1,hc=0.7),
                w=c(0,10,0,10))
   X3.straush <- rmh(model=mod06,start=list(n.start=ns),
                     control=list(nrep=nr,nverb=nv))
   
   # Soft core:
   w    <- c(0,10,0,10)
   mod07 <- list(cif="sftcr",par=list(beta=0.8,sigma=0.1,kappa=0.5),
                w=c(0,10,0,10))
   X.sftcr <- rmh(model=mod07,start=list(n.start=ns),
                  control=list(nrep=nr,nverb=nv))
   if(interactive()) plot(X.sftcr)

   # Area-interaction process:
   mod42 <- rmhmodel(cif="areaint",par=list(beta=2,eta=1.6,r=0.7),
                 w=c(0,10,0,10))
   X.area <- rmh(model=mod42,start=list(n.start=ns),
                  control=list(nrep=nr,nverb=nv))
   if(interactive()) plot(X.area)

   # Triplets process
   modtrip <- list(cif="triplets",par=list(beta=2,gamma=0.2,r=0.7),
                   w=c(0,10,0,10))
   X.triplets <- rmh(model=modtrip,
                     start=list(n.start=ns),
                     control=list(nrep=nr,nverb=nv))
   if(interactive()) plot(X.triplets)
   
   # Multitype Strauss:
   beta <- c(0.027,0.008)
   gmma <- matrix(c(0.43,0.98,0.98,0.36),2,2)
   r    <- matrix(c(45,45,45,45),2,2)
   mod08 <- list(cif="straussm",par=list(beta=beta,gamma=gmma,radii=r),
                w=c(0,250,0,250))
   X1.straussm <- rmh(model=mod08,start=list(n.start=ns),
                      control=list(ptypes=c(0.75,0.25),nrep=nr,nverb=nv))
   if(interactive()) plot(X1.straussm)
   
   # Multitype Strauss conditioning upon the total number
   # of points being 80:
   X2.straussm <- rmh(model=mod08,start=list(n.start=ns),
                      control=list(p=1,ptypes=c(0.75,0.25),nrep=nr,
                                   nverb=nv))
   
   # Conditioning upon the number of points of type 1 being 60
   # and the number of points of type 2 being 20:
   X3.straussm <- rmh(model=mod08,start=list(n.start=c(60,20)),
                      control=list(fixall=TRUE,p=1,ptypes=c(0.75,0.25),
                                   nrep=nr,nverb=nv))
   
   # Multitype Strauss hardcore:
   rhc  <- matrix(c(9.1,5.0,5.0,2.5),2,2)
   mod09 <- list(cif="straushm",par=list(beta=beta,gamma=gmma,
                iradii=r,hradii=rhc),w=c(0,250,0,250))
   X.straushm <- rmh(model=mod09,start=list(n.start=ns),
                     control=list(ptypes=c(0.75,0.25),nrep=nr,nverb=nv))
   
   # Multitype Strauss hardcore with trends for each type:
   beta  <- c(0.27,0.08)
   tr3   <- function(x,y){x <- x/250; y <- y/250;
   			   exp((6*x + 5*y - 18*x^2 + 12*x*y - 9*y^2)/6)
                         }
                         # log quadratic trend
   tr4   <- function(x,y){x <- x/250; y <- y/250;
                         exp(-0.6*x+0.5*y)}
                        # log linear trend
   mod10 <- list(cif="straushm",par=list(beta=beta,gamma=gmma,
                 iradii=r,hradii=rhc),w=c(0,250,0,250),
                 trend=list(tr3,tr4))
   X1.straushm.trend <- rmh(model=mod10,start=list(n.start=ns),
                            control=list(ptypes=c(0.75,0.25),
                            nrep=nr,nverb=nv))
   if(interactive()) plot(X1.straushm.trend)
   
   # Multitype Strauss hardcore with trends for each type, given as images:
   bigwin <- square(250)
   i1 <- as.im(tr3, bigwin)
   i2 <- as.im(tr4, bigwin)   
   mod11 <- list(cif="straushm",par=list(beta=beta,gamma=gmma,
                 iradii=r,hradii=rhc),w=bigwin,
                 trend=list(i1,i2))
   X2.straushm.trend <- rmh(model=mod11,start=list(n.start=ns),
                            control=list(ptypes=c(0.75,0.25),expand=1,
                            nrep=nr,nverb=nv))
   
   # Diggle, Gates, and Stibbard:
   mod12 <- list(cif="dgs",par=list(beta=3600,rho=0.08),w=c(0,1,0,1))
   X.dgs <- rmh(model=mod12,start=list(n.start=ns),
                control=list(nrep=nr,nverb=nv))
   if(interactive()) plot(X.dgs)
   
   # Diggle-Gratton:
   mod13 <- list(cif="diggra",
                 par=list(beta=1800,kappa=3,delta=0.02,rho=0.04),
                 w=square(1))
   X.diggra <- rmh(model=mod13,start=list(n.start=ns),
                   control=list(nrep=nr,nverb=nv))
   if(interactive()) plot(X.diggra)
   
   # Fiksel:
   modFik <- list(cif="fiksel",
                 par=list(beta=180,r=0.15,hc=0.07,kappa=2,a= -1.0),
                 w=square(1))
   X.fiksel <- rmh(model=modFik,start=list(n.start=ns),
                   control=list(nrep=nr,nverb=nv))
   if(interactive()) plot(X.fiksel)
   
   # Geyer:
   mod14 <- list(cif="geyer",par=list(beta=1.25,gamma=1.6,r=0.2,sat=4.5),
                 w=c(0,10,0,10))
   X1.geyer <- rmh(model=mod14,start=list(n.start=ns),
                   control=list(nrep=nr,nverb=nv))
   if(interactive()) plot(X1.geyer)
   
   # Geyer; same as a Strauss process with parameters
   # (beta=2.25,gamma=0.16,r=0.7):
   
   mod15 <- list(cif="geyer",par=list(beta=2.25,gamma=0.4,r=0.7,sat=10000),
                 w=c(0,10,0,10))
   X2.geyer <- rmh(model=mod15,start=list(n.start=ns),
                   control=list(nrep=nr,nverb=nv))
   
   mod16 <- list(cif="geyer",par=list(beta=8.1,gamma=2.2,r=0.08,sat=3))

   X3.geyer <- rmh(model=mod16,start=list(x.start=redwood),
                   control=list(periodic=TRUE,nrep=nr,nverb=nv))
   
   # Geyer, starting from the redwood data set, simulating
   # on a torus, and conditioning on n:
   X4.geyer <- rmh(model=mod16,start=list(x.start=redwood),
                   control=list(p=1,periodic=TRUE,nrep=nr,nverb=nv))

   # Lookup (interaction function h_2 from page 76, Diggle (2003)):
      r <- seq(from=0,to=0.2,length=101)[-1] # Drop 0.
      h <- 20*(r-0.05)
      h[r<0.05] <- 0
      h[r>0.10] <- 1
      mod17 <- list(cif="lookup",par=list(beta=4000,h=h,r=r),w=c(0,1,0,1))
      X.lookup <- rmh(model=mod17,start=list(n.start=ns),
                      control=list(nrep=nr,nverb=nv))
      if(interactive()) plot(X.lookup)
                   
   # Strauss with trend
   tr <- function(x,y){x <- x/250; y <- y/250;
   			   exp((6*x + 5*y - 18*x^2 + 12*x*y - 9*y^2)/6)
                         }
   beta <- 0.3
   gmma <- 0.5
   r    <- 45
   modStr <- list(cif="strauss",par=list(beta=beta,gamma=gmma,r=r),
                 w=square(250), trend=tr)
   X1.strauss.trend <- rmh(model=modStr,start=list(n.start=ns),
                           control=list(nrep=nr,nverb=nv))
   # Baddeley-Geyer
   r <- seq(0,0.2,length=8)[-1]
   gmma <- c(0.5,0.6,0.7,0.8,0.7,0.6,0.5)
   mod18 <- list(cif="badgey",par=list(beta=4000, gamma=gmma,r=r,sat=5),
                 w=square(1))
   X1.badgey <- rmh(model=mod18,start=list(n.start=ns),
                    control=list(nrep=nr,nverb=nv))
   mod19 <- list(cif="badgey",
                 par=list(beta=4000, gamma=gmma,r=r,sat=1e4),
                 w=square(1))
   set.seed(1329)
   X2.badgey <- rmh(model=mod18,start=list(n.start=ns),
                    control=list(nrep=nr,nverb=nv))

   # Check:
   h <- ((prod(gmma)/cumprod(c(1,gmma)))[-8])^2
   hs <- stepfun(r,c(h,1))
   mod20 <- list(cif="lookup",par=list(beta=4000,h=hs),w=square(1))
   set.seed(1329)
   X.check <- rmh(model=mod20,start=list(n.start=ns),
                      control=list(nrep=nr,nverb=nv))
   # X2.badgey and X.check will be identical.

   mod21 <- list(cif="badgey",par=list(beta=300,gamma=c(1,0.4,1),
                 r=c(0.035,0.07,0.14),sat=5), w=square(1))
   X3.badgey <- rmh(model=mod21,start=list(n.start=ns),
                    control=list(nrep=nr,nverb=nv))
   # Same result as Geyer model with beta=300, gamma=0.4, r=0.07,
   # sat = 5 (if seeds and control parameters are the same)

   # Or more simply:
   mod22 <- list(cif="badgey",
                 par=list(beta=300,gamma=0.4,r=0.07, sat=5),
                 w=square(1))
   X4.badgey <- rmh(model=mod22,start=list(n.start=ns),
                    control=list(nrep=nr,nverb=nv))
   # Same again --- i.e. the BadGey model includes the Geyer model.


   # Illustrating scalability.
   if(FALSE) {
    M1 <- rmhmodel(cif="strauss",par=list(beta=60,gamma=0.5,r=0.04),w=owin())
    set.seed(496)
    X1 <- rmh(model=M1,start=list(n.start=300))
    M2 <- rmhmodel(cif="strauss",par=list(beta=0.6,gamma=0.5,r=0.4),
              w=owin(c(0,10),c(0,10)))
    set.seed(496)
    X2  <- rmh(model=M2,start=list(n.start=300))
    chk <- affine(X1,mat=diag(c(10,10)))
    all.equal(chk,X2,check.attributes=FALSE)
    # Under the default spatstat options the foregoing all.equal()
    # will yield TRUE.  Setting spatstat.options(scalable=FALSE) and
    # re-running the code will reveal differences between X1 and X2.
   }

   if(!interactive()) spatstat.options(oldopt)

---