spatstat (version 1.3-4)

rmh.default: 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: [object Object],[object Object],[object Object],[object Object] See Details for details.
start
List of parameters determining the initial state of the algorithm: [object Object],[object Object],[object Object] 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: [object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object] See Details for d
...
Further arguments; currently ignored.

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.

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

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:

[object Object],[object Object],[object Object],[object Object] The argument start determines the initial state of the Metropolis-Hastings algorithm. Possible components are [object Object],[object Object],[object Object] 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: [object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object]

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{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. }

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

  • The value ofexpandmust be equal to 1; it defaults to 1, and it is an error to specify a value larger than 1.
The window in which the pattern is being simulated must (currently) be rectangular. The resulting simulated pattern will have precisely the number of points (of each type) specified by n.start. }

See Also

rmh, rmh.ppm, ppp, mpl, Strauss, Softcore, StraussHard, MultiStrauss, MultiStraussHard, DiggleGratton

Examples

Run this code
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))
   <testonly># 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=100,nverb=50))

    # Strauss process, conditioning on n = 80:
    X2.strauss <- rmh(model=mod01,start=list(n.start=80),
                      control=list(p=1,nrep=10,nverb=5))
    
    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=10,nverb=5))
    

    # Strauss process == 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=10,nverb=5,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=10,nverb=5))
    
    # 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=10,nverb=5))
    
    # 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=10,nverb=5))
    
    # 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=10,nverb=5,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=10,nverb=5))
    
    # 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=10,nverb=5))
    
    # 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=10,nverb=5))
    
    # 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=10,nverb=5))

    # 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=10,nverb=5))
    
    # 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=10,nverb=5))
    
    # 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=10,nverb=5))
    
    # 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=10,nverb=5))
    
    # 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=10,nverb=5))

    # 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=10,nverb=5))

    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=10,nverb=5))

    # 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=10,nverb=5))</testonly>

Run the code above in your browser using DataCamp Workspace