rmh.default

0th

Percentile

Simulate Point Process Models using the Metropolis-Hastings Algorithm.

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

Keywords
spatial, datagen
Usage
## S3 method for class 'default':
rmh(model, start, control, verbose=TRUE,
track=FALSE, \dots)
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.
verbose
Logical flag indicating whether to print progress reports.
track
Logical flag indicating whether to save the transition history of the simulations.
...
Further arguments which will be passed to any trend functions in model.
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:

[object Object],[object Object],[object Object],[object Object],[object Object] For full details of these parameters, see rmhmodel.

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

[object Object],[object Object] 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 a list, or an object of class "rmhcontrol", with components: [object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object] For full details of these parameters, see rmhcontrol.

Value

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

    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 value of .Random.seed at the start of the simulations is also saved and returned as an attribute seed.

    If track=TRUE, 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.

Conditional Simulation

There are several kinds of conditional simulation.

  • Simulationconditional upon the number of points, that is, holding the number of points fixed. To do this, setcontrol$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 settingstart$n.startto be a scalar, or by setting the initial patternstart$x.start.
  • In the case of multitype processes, it is possible to simulate the modelconditionally upon the number of points of each type, i.e. holding the number of points of each type to be fixed. To do this, setcontrol$pequal to 1 andcontrol$fixallto beTRUE. The number of points is then determined by the starting state, which may be specified either by settingstart$n.startto be an integer vector, or by setting the initial patternstart$x.start.
  • Simulationconditional 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, setcontrol$x.condto equal the specified point pattern$y$, making sure that it is an object of class"ppp"and that the windowas.owin(control$x.cond)is the conditioning window$V$.
  • Simulationconditional 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, setcontrol$x.condto be adata.framecontaining 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.

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 x.start$window. Any specified value for expand is simply ignored.

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

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.

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, AreaInter, BadGey, DiggleGatesStibbard, DiggleGratton, Fiksel, Geyer, Hardcore, LennardJones, MultiHard, MultiStrauss, MultiStraussHard, PairPiece, Poisson, Softcore, Strauss, StraussHard

Aliases
  • rmh.default
Examples
nr   <- 1e5
   nv  <- 5000
   ns <- 200
   <testonly>nr  <- 10
   nv <- 5
   ns <- 20
   oldopt <- spatstat.options()
   spatstat.options(expand=1.1)</testonly>
   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))
   
   # Strauss process, conditioning on n = 80:
   X2.strauss <- rmh(model=mod01,start=list(n.start=ns),
                     control=list(p=1,nrep=nr,nverb=nv))

   <testonly>X1.strauss <- rmh(model=mod01,start=list(n.start=ns),
                     control=list(nrep=nr,nverb=nv), track=TRUE)
     X2.strauss <- rmh(model=mod01,start=list(n.start=ns),
                     control=list(p=1,nrep=nr,nverb=nv), track=TRUE)</testonly>
   
   # 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))
   
   # 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))
   
   # 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
   xxx$window <- 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))

   # 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))

   
   # 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))
   
   # 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))
   
   # 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))
   
   # 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))
   
   # Fiksel:
   modFik <- list(cif="fiksel",
                 par=list(beta=180,r=0.15,hc=0.07,kappa=2,a= -1.0),
                 w=square(1))
   X.diggra <- rmh(model=modFik,start=list(n.start=ns),
                   control=list(nrep=nr,nverb=nv))
   
   # 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))
   
   # 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))
   data(redwood)
   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))
                   
   # 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
   mod17 <- list(cif="strauss",par=list(beta=beta,gamma=gmma,r=r),
                 w=square(1), trend=tr3)
   X1.strauss.trend <- rmh(model=mod17,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.
   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.attribute=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.

   <testonly>spatstat.options(oldopt)</testonly>
Documentation reproduced from package spatstat, version 1.23-3, License: GPL (>= 2)

Community examples

Looks like there are no examples yet.