rmh.default {spatstat}R Documentation

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

rmh.default(model,start,control,...) 

Arguments

model List of parameters specifying the point process model that is to be simulated:
cif
Character string specifying the choice of model
par
Parameters of the model
tpar
(optional) Parameters of the spatial trend
w
Spatial window in which to simulate
See Details for details.
start List of parameters determining the initial state of the algorithm:
n.start
number of initial points (randomly generated)
x.start
initial point pattern configuration
iseed
vector of seeds for random number generator
The parameters n.start and x.start are incompatible. See Details for details.
control List of parameters controlling the iterative behaviour and termination of the algorithm:
expand
Factor by which the window will be expanded
periodic
(Logical) whether to apply periodic boundary conditions
nrep
Total number of steps (proposals) of Metropolis-Hastings algorithm that should be run
p
Probability of proposing a shift (as against a birth/death)
q
Conditional probability of proposing a death given that a birth or death will be proposed
ptypes
For multitype point processes, the distribution of the mark attached to a new random point (when a birth is proposed)
fixall
(Logical) for multitype point processes, whether to fix the number of points of each type.
nverb
Progress reports will be printed every nverb iterations
See Details for details.
... Further arguments; currently ignored.

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 a list with the following components:

cif
A character string specifying the choice of interpoint interaction for the point process. The current options are
'strauss'
The Strauss process
'straush'
The Strauss process with hard core
'sftcr'
The Softcore process
'straussm'
The multitype Strauss process
'straushm'
Multitype Strauss process with hard core
'dgs'
Diggle, Gates and Stibbard (1987) process
'diggra'
Diggle and Gratton (1984) process
'geyer'
Saturation process (Geyer, 1999).

See the section Extensions for the possibility of extending this list of options.

par
A vector (or a list if the pattern is multitype) providing a set of parameter values appropriate to the conditional intensity function being invoked. par should be as follows, for each of the available conditional intensity functions:
strauss:
(Strauss process.) A named vector with components beta,gamma,r which are respectively the ``base'' intensity, the pair-wise interaction parameter and the interaction radius. Note that gamma must be less than or equal to 1.
straush:
(Strauss process with hardcore.) A named vector with entries beta,gamma,r,hc where beta, gamma, and r are as for the Strauss process, and hc is the hardcore radius. Of course hc must be less than r.
sftcr:
(Softcore process.) A named vector with components beta,sigma,kappa. Again beta is a ``base'' intensity. The pairwise interaction between two points u != v is

-(sigma/||u-v||)^(2/kappa)

Note that it is necessary that 0 < kappa <1.

straussm:
(Multitype Strauss process.) A named list with components
  • beta: A vector of ``base'' intensities, one for each possible type.
  • gamma: A symmetric matrix of interaction parameters, with gamma_ij pertaining to the interaction between type i and type j.
  • radii: A symmetric matrix of interaction radii, with entries r_ij pertaining to the interaction between type i and type j.
straushm:
(Multitype Strauss process with hardcore.) A named list with components beta and gamma as for straussm and two ``radii'' components:
  • iradii: the interaction radii
  • hradii: the hardcore radii

which are both symmetric matrices of nonnegative numbers. The entries of hradii must be less than the corresponding entries of iradii.

dgs:
(Diggle, Gates, and Stibbard process. See Diggle, Gates, and Stibbard (1987)) A named vector with components beta and rho. This process has pairwise interaction function

e(t) = sin^2((pi * t)/(2 * rho))

for t < rho, and equal to 1 for t >= rho.

diggra:
(Diggle-Gratton process. See Diggle and Gratton (1984) and Diggle, Gates and Stibbard (1987).) A named vector with components beta, kappa, delta and rho. This process has pairwise interaction function e(t) equal to 0 for t < delta, equal to

((t-delta)/(rho-delta))^kappa

for delta <= t < rho, and equal to 1 for delta >= rho. Note that here we use the symbol kappa where Diggle, Gates, and Stibbard use beta since we reserve the symbol beta for an intensity parameter.

geyer
(See Geyer (1999)) A named vector with components beta, gamma, r, and sat. The components beta, gamma, r are as for the Strauss model, and sat is the ``saturation'' parameter. The model is Geyer's ``saturation'' point process model. It is ``like a Strauss model, but with an upper bound to the number of r-close neighbors of any point.''

Explicitly, a saturation point process with interaction radius r, saturation threshold s, and parameters beta and gamma, is the point process in which each point x[i] in the pattern X contributes a factor

beta gamma^min(s,t(x[i],X))

to the probability density of the point pattern, where t(x[i],X) denotes the number of ``r-close neighbours'' of x[i] in the pattern X.

tpar
A vector (or a list if the pattern is multitype) specifying the coefficients of a polynomial for a log polynomial trend. The coefficients (corresponding to each type, for a multitype process) must be given in the order of the terms x, y, x^2, xy, y^2, x^3, x^2y, xy^2, y^3, ..., y^n where n is the degree of the trend. (Thus the possible lengths of a sequence of coefficients are 0, 2, 5, 9, 14, ...). The degree of the polynomial will be calculated from the length of tpar. If the user wishes to omit some of the terms in such a polynomial, the corresponding coefficients should be specified as 0. For example, to specify a trend 2x the trend parmeters should be specified as c(2,0).

If the process is unmarked tpar must be a vector consisting solely of the coefficients of the trend polynomial in the correct order.

If the process is multitype, so that tpar is given as a list, the i-th component of this list must consist of the coefficients for the trend corresponding to the i-th type. If there is to be no trend for the i-th type, then the i-th component of the list must be NULL (but must be present in the list as the NULL object.)

w
A specification of a window in which the pattern is to be generated. If specified, it must be in a form which can be coerced to an object of class owin by as.owin(). Note that for non-rectangular windows we cannot (as yet) generate a pattern from a process which exists only within that window. The domain of the theoretical process must (at present) be at least as large as the ``enclosing box'' of the window.

The argument start determines the initial state of the Metropolis-Hastings algorithm. Possible components are

n.start
The number of ``initial'' points to be randomly (uniformly) generated in the window w. For a multitype point process, n.start may be a vector (of length equal to the number of types) giving the number of points of each type to be generated. Incompatible with x.start.

n.start is the number of ``initial'' points to be randomly (uniformly) generated in the window w, or a vector (of length equal to the number of types for a multitype process) giving the number of points of each type to be generated. These uniformly generated points form an initial state or configuration for the Metropolis-Hastings algorithm. Precisely one of n.start or x.start must be specified.

A vector-valued n.start is meaningful only if p (the probability of a shift as opposed to a birth or death) is equal to 1 (so that we are conditioning on the number of points). When p < 1, if n.start is vector valued then it is effectively replaced by its sum.

The resulting set of uniformly generated points gives the Metropolis-Hastings algorithm an initial state from which to start. (Actually, when p < 1, the number n.start gets multiplied by the ratio of the area of the enclosing box for the window to the area of the window, and then by the factor expand. Then that many points are uniformly generated in the expanded window; see below.) The value of n.start should be roughly equal to (an educated guess at) the expected number of points which will be generated inside the window.

x.start
Initial point pattern configuration. A point pattern (an object of class "ppp", or data which can be coerced to this class by as.ppp). Incompatible with n.start.

x.start is a point pattern (an object of class ppp, or an object which can be coerced to this class by as.ppp()). This object provides an alternative means of specifying the initial ``state'' or configuration for the Metropolis-Hastings algorithm. If x.start is specified, but x.start$window is NULL, then this gap is filled in by the component w of start. If x.start$window is present, and if w is specified as well, the latter is used as a window to which to clip the final simulated pattern. Thus in such circumstances it is only sensible to specify a value of w which is contained in x.start$window. However no checking is done for this.

The simulated pattern is constructed in the enclosing rectangle of x.start$window. No expansion takes place. (The argument expand is forced to equal 1.) As indicated above, at the end of the simulation, the resulting pattern is clipped to the window w if this is given (or to x.start$window if this differs from its enclosing rectangle).

iseed
Seed for random number generator. A triple of integers. If unspecified these are themselves generated, on the interval from 1 to 1 million, using the function sample.

The parameters n.start and x.start are incompatible.

The third argument control controls the simulation procedure, iterative behaviour, and termination of the Metropolis-Hastings algorithm. It is a list with components:

expand
The factor by which the enclosing box of the window w is to be expanded in order to better approximate the simulation of a process existing in the whole plane, rather than just in the enclosing box. If expand equals 1, then we are simulating the latter (unless periodic is TRUE; see below). The larger expand is, the better we approximate the former. Note that any value of expand smaller than 1 is treated as if it were 1.

The area of the expanded window is equal to expand times the area of the enclosing box; width and height are stretched proportionately. Points are generated by the Metropolis-Hastings algorithm in the expanded window, and then ``clipped'' down to the original window w (or to x.start$window) when the algorithm has finished. The argument expand defaults to 2 if periodic is FALSE and p < 1 and to 1 if periodic is TRUE or if p < 1, or if the starting configuration is specified via x.start. Trying to set expand greater than 1 when periodic is TRUE or p = 1 generates an error. A specified value of expand is simply ignored if x.start is given.

periodic
A logical scalar; if periodic is TRUE we simulate a process on the torus formed by identifying opposite edges of the (rectangular) window. If periodic is TRUE and the window (i.e. x.start$window if x.start is specified; otherwise w) is not rectangular, an error is given.
nrep
The number of repetitions or steps (changes of state) to be made by the Metropolis-Hastings algorithm. It should be large.
p
The probability of of proposing a ``shift'' (as opposed to a birth or death) in the Metropolis-Hastings algorithm. If p = 1 then we do nothing but shifts, whence the number of points never changes, i.e. we are simulating conditionally upon the number of points. In this case, for multitype processes, we also have the option of simulating conditionally upon the number of points of each type; this is effected by specifying fixall=TRUE. In this case, n.start must be a vector whose entries are these numbers.

We can only condition on the number of points if the simulation takes place in the original window (as opposed to taking place in a larger window and then being clipped to the original). Hence, if p = 1 then expand defaults to 1 and it is an error to specify a value of expand which is greater than 1. Likewise, if p = 1, it is (currently) an error to specify anything other than a rectangular window.

q
The probability of proposing a death (rather than a birth) given that birth/death has been chosen over shift. This is of course ignored if p is equal to 1.
ptypes
A vector of probabilities (summing to 1) to be used in assigning a random type to a new point. Defaults to a vector each of whose entries is 1/nt where nt is the number of types for the process. Convergence of the simulation algorithm should be improved if ptypes is close to the relative frequencies of the types which will result from the simulation.
fixall
A logical scalar specifying whether to condition on the number of points of each type. Meaningful only if a marked process is being simulated, and if p = 1. A warning message is given if fixall is set equal to TRUE when it is not meaningful.
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.

Value

A point pattern (an object of class "ppp", see ppp.object).
The return value has an attribute info consisting of arguments supplied to the function (or default values of arguments which were not explicitly supplied). These are given so that it is possible to reconstruct exactly the manner in which the pattern was generated. The components of info are model, start, and control which in turn are lists:
model=list(cif, par, tpar)
start=list(n.start,x.start,iseed)
control=list(expand,periodic,nrep,p,q,ptypes,fixall,nverb).
Note that if x.start is specified only its name is preserved inside info.

Note

It is possible to simulate conditionally upon the number of points, or in the case of multitype processes, upon the number of points of each type. To condition upon the total number of points, set p (the probability of a shift) equal to 1, and specify n.start to be a scalar (as usual). To condition upon the number of points of each type, set p equal to 1, fixall equal to TRUE, and specify n.start to be a vector of length nt where nt is the number of types.

In these circumstances

Warnings

There is never a guarantee that the Metropolis-Hastings algorithm has converged to the steady state.

If x.start is specified then expand is set equal to 1 and simulation takes place in the enclosing rectangle of x.start$window. Any specified value for expand is simply ignored.

If trends are specified, make sure that the lengths of the vectors of coefficients in tpar make sense. 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.

Extensions

The argument model$cif matches the name of a Fortran subroutine which calculates the conditional intensity function for the model. It is intended that more options will be added in the future. The very brave user could try to add her own. Note that in addition to writing Fortran code for the new conditional intensity function, the user would have to modify the code in the files cif.f and rmh.default.R appropriately. (And of course re-install the spatstat package so as to update the dynamically loadable shared object spatstat.so.)

Author(s)

Adrian Baddeley adrian@maths.uwa.edu.au http://www.maths.uwa.edu.au/~adrian/ and Rolf Turner rolf@math.unb.ca http://www.math.unb.ca/~rolf

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. 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 M{o}ller, 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, ppp, mpl, Strauss, Softcore, StraussHard, MultiStrauss, MultiStraussHard, DiggleGratton

Examples

   library(spatstat)
   
    # Strauss process.
    mod01 <- list(cif="strauss",par=c(beta=2,gamma=0.2,r=0.7),
                  w=c(0,10,0,10))
    X1.strauss <- rmh(model=mod01,start=list(n.start=80),
                      control=list(nrep=1e5,nverb=5000))

    # Strauss process, conditioning on n = 80:
    X2.strauss <- rmh(model=mod01,start=list(n.start=80),
                      control=list(p=1,nrep=1e5,nverb=5000))
    
    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)
    mod02 <- list(cif="strauss",par=c(beta=2000,gamma=0.6,r=0.7),
                 w=owin(poly=list(x=x,y=y)))
    X3.strauss <- rmh(model=mod02,start=list(n.start=90),
                      control=list(nrep=1e5,nverb=5000))

    # Strauss process equal to pure hardcore:
    mod03 <- list(cif="strauss",par=c(beta=2,gamma=0,r=0.7),w=c(0,10,0,10))
    X4.strauss <- rmh(model=mod03,start=list(n.start=60),
                      control=list(nrep=1e5,nverb=5000,iseed=c(42,17,69)))

    # Another Strauss process, starting off from X3.strauss (but
    # with a rectangular window):
    xxx <- X3.strauss
    xxx$w <- as.owin(c(0,1,0,1))
    X3.strauss <- rmh(model=mod02,start=list(x.start=xxx),
                      control=list(nrep=1e5,nverb=5000))
    
    # Strauss with hardcore:
    mod04 <- list(cif="straush",par=c(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=70),
                      control=list(nrep=1e5,nverb=5000))
    
    # Another Strauss with hardcore (with a perhaps surprizing
    # result):
    mod05 <- list(cif="straush",par=c(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=250),
                      control=list(nrep=1e5,nverb=5000))
    
    # Pure hardcore (identical to X3.strauss).
    mod06 <- list(cif="straush",par=c(beta=2,gamma=1,r=1,hc=0.7),
                 w=c(0,10,0,10))
    X3.straush <- rmh(model=mod06,start=list(n.start=60),
                      control=list(nrep=1e5,nverb=5000,iseed=c(42,17,69)))
    
    # Soft core:
    par3 <- c(0.8,0.1,0.5)
    w    <- c(0,10,0,10)
    mod07 <- list(cif="sftcr",par=c(beta=0.8,sigma=0.1,kappa=0.5),
                 w=c(0,10,0,10))
    X.sftcr <- rmh(model=mod07,start=list(n.start=70),
                   control=list(nrep=1e5,nverb=5000))
    
    # 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=80),
                       control=list(ptypes=c(0.75,0.25),nrep=1e5,nverb=5000))
    
    # Multitype Strauss conditioning upon the total number
    # of points being 80:
    X2.straussm <- rmh(model=mod08,start=list(n.start=80),
                       control=list(p=1,ptypes=c(0.75,0.25),nrep=1e5,
                                    nverb=5000))
    
    # 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=1e5,nverb=5000))

    # 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=80),
                      control=list(ptypes=c(0.75,0.25),nrep=1e5,nverb=5000))
    
    # Multitype Strauss hardcore with trends for each type:
    beta  <- c(0.0027,0.08)
    tpar1 <- c(0.02,0.004,-0.0004,0.004,-0.0004) # Coefs. for log quadratic
    tpar2 <- c(-0.06,0.05)                       # and log linear trends.
    mod10 <- list(cif="straushm",par=list(beta=beta,gamma=gmma,
                  iradii=r,hradii=rhc),w=c(0,250,0,250),
                  tpar=list(tpar1,tpar2))
    X1.straushm.trend <- rmh(model=mod10,start=list(n.start=350),
                             control=list(ptypes=c(0.75,0.25),
                             nrep=1e5,nverb=5000))
    
    # Diggle, Gates, and Stibbard:
    mod11 <- list(cif="dgs",par=c(beta=3600,rho=0.08),w=c(0,1,0,1))
    X.dgs <- rmh(model=mod11,start=list(n.start=300),
                 control=list(nrep=1e5,nverb=5000))
    
    # Diggle-Gratton:
    mod12 <- list(cif="diggra",
                  par=c(beta=1800,kappa=3,delta=0.02,rho=0.04),
                  w=square(1))
    X.diggra <- rmh(model=mod12,start=list(n.start=300),
                    control=list(nrep=1e5,nverb=5000))
    
    # Geyer:
    mod13 <- list(cif="geyer",par=c(beta=1.25,gamma=1.6,r=0.2,sat=4.5),
                  w=c(0,10,0,10))
    X1.geyer <- rmh(model=mod13,start=list(n.start=200),
                    control=list(nrep=1e5,nverb=5000))

    # Geyer; same as a Strauss process with parameters
    # (beta=2.25,gamma=0.16,r=0.7):

    mod14 <- list(cif="geyer",par=c(beta=2.25,gamma=0.4,r=0.7,sat=10000),
                  w=c(0,10,0,10))
    X2.geyer <- rmh(model=mod14,start=list(n.start=200),
                    control=list(nrep=1e5,nverb=5000))

    mod15 <- list(cif="geyer",par=c(beta=8.1,gamma=2.2,r=0.08,sat=3))
    data(redwood)
    X3.geyer <- rmh(model=mod15,start=list(x.start=redwood),
                    control=list(periodic=TRUE,nrep=1e5,nverb=5000))

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

[Package Contents]