rmh: Simulate point patterns using the Metropolis-Hastings algorithm.

Description

Simulates point patterns from a (still small) range of point process models.

Usage

rmh(cif,par,w=NULL,ntypes=1,ptypes=NULL,tpar=NULL,n.start=NULL,
        x.start=NULL,fixall=FALSE,expand=NULL,periodic=F,nrep=1e6,
        p=0.9,q=0.5,iseed=NULL,nverb=0)

Arguments

cif

A character string giving the name of a Fortran subroutine which will calculate the value of the conditional intensity function for the model. So far the only permissible values are strauss (Strauss process), str

par

A set of parameter values appropriate to the conditional intensity function being invoked. (See Details.)

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

ntypes

The number of distinct types to be used if the simulated pattern is to be a multitype (marked) point pattern.

ptypes

A vector of probabilities (of length ntypes and 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/ntypes. Convergence of the simulation algorithm sho

tpar

A vector (or possibly 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

n.start

The number of ``initial'' points to be randomly (uniformly) generated in the window w, or a vector (of length ntypes) giving the number of points of each type to be generated. These uniformly generated points form an

x.start

A point pattern (an object of class ppp), which provides an alternative means of specifying the initial ``state'' or configuration for the Metropolis-Hastings algorithm. The simulated pattern is constructed in the enclosing rectan

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 w

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

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 (x.start$w

nrep

The number of repetitions or steps (changes of state) to be made by the Metropolis-Hastings algorithm. It should be large.

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

iseed

A vector equal to a triple of integers to be used as seeds to the random number generating procedure. If unspecified these are themselves generated, on the interval from 1 to 1 million, using the function sample().

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 list of class ppp with the usual components
windowThe window w in the form of an object of class owin.
nThe number of generated points (in the window w).
xThe x-coordinates of the generated points.
yThe y-coordinates of the generated points.
marksThe marks of the generated points. Remember that this component is a factor. (This component is present only when the model is a multitype point process.)
In addition the list returned has a component 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: cif, par, tpar, n.start, x.start.name, nrep, p, q, expand, periodic, and iseed.
Note that if x.start is specified only its name is preserved in the list 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$w. 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, it is probably safest to specify tpar as a list. But 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.

Note that if tpar is given as an atomic vector, for multitype processes, no checking is (or, realistically, can be) done to make sure that the degrees supplied are sensible.

Details

The parameter vector (or list) par should be as follows, for each of the available conditional intensity functions: [object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object]

References

Baddeley, Adrian, and Turner, Rolf. ``Practical maximum pseudolikelihood for spatial point patterns.'' Australian and New Zealand Journal of Statistics, vol. 42, 2000, pp. 283 -- 322.

Diggle, Peter J., Gates, David J., and Stibbard, Alyson. ``A nonparametric estimator for pairwise-interaction point processes.'' Biometrika, vol. 74, 1987, pp. 763 -- 770.

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

Examples

Run this code

library(spatstat)
   par11 <- c(2,0.2,0.7)
    w     <- c(0,10,0,10)
    X1.strauss <- rmh("strauss",par=par11,w=w,n.start=80,
                      nrep=1e5,nverb=5000)

    # Conditioning on n = 80:
    X2.strauss <- rmh("strauss",par=par11,w=w,n.start=80,
                      p=1,nrep=1e5,nverb=5000)
    
    par12 <- c(2000,0.6,0.07)
    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)
    w     <- owin(poly=list(x=x,y=y))
    X3.strauss <- rmh("strauss",par=par12,w=w,n.start=90,
                      nrep=1e5,nverb=5000)
    
    # Pure hardcore:
    par13 <- c(2,0,0.7)
    w     <- c(0,10,0,10)
    X4.strauss <- rmh("strauss",par=par13,w=w,n.start=60,
                      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))
    X5.strauss <- rmh("strauss",par=par12,w=w,x.start=xxx,
                      nrep=1e5,nverb=5000)
    
    par21 <- c(2,0.2,0.7,0.3)
    w     <- c(0,10,0,10)
    X1.straush <- rmh("straush",par=par21,w=w,n.start=70,
                      nrep=1e5,nverb=5000)
    
    par22 <- c(80,0.36,45,2.5)
    w     <- c(0,250,0,250)
    X2.straush <- rmh("straush",par=par22,w=w,n.start=160,
                      nrep=1e5,nverb=5000)
    
    # Pure hardcore (identical to X3.strauss).
    par23 <- c(2,1,1,0.7)
    w     <- c(0,10,0,10)
    X3.straush <- rmh("straush",par=par23,w=w,n.start=60,
                      nrep=1e5,nverb=5000,iseed=c(42,17,69))
    
    par3 <- c(0.8,0.1,0.5)
    w    <- c(0,10,0,10)
    X.sftcr <- rmh("sftcr",par=par3,w=w,n.start=70,nrep=1e5,
                   nverb=5000)
    
    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)
    par4 <- list(beta=beta,gamma=gmma,r=r)
    w    <- c(0,250,0,250)
    pm   <- c(0.75,0.25)
    X1.straussm <- rmh("straussm",par=par4,w=w,ntypes=2,
                      ptypes=pm,n.start=80,nrep=1e5,nverb=5000)
    
    # Conditioning upon the total number of points being 80:
    X2.straussm <- rmh("straussm",par=par4,w=w,ntypes=2,
                      ptypes=pm,n.start=80,p=1,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("straussm",par=par4,w=w,ntypes=2,n.start=c(60,20),
                       fixall=TRUE,nrep=1e5,p=1,nverb=5000)

    rhc  <- matrix(c(9.1,5.0,5.0,2.5),2,2)
    par5 <- list(beta=beta,gamma=gmma,r=r,rhc=rhc)
    X.straushm <- rmh("straushm",par=par5,w=w,ntypes=2,
                      ptypes=pm,n.start=80,nrep=1e5,nverb=5000)
    
    beta  <- c(0.0027,0.08)
    par6  <- list(beta=beta,gamma=gmma,r=r,rhc=rhc)
    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.
    w     <- c(0,250,0,250)
    pm    <- c(0.75,0.25)
    X1.straushm.trend <- rmh("straushm",par=par6,w=w,ntypes=2,
                             ptypes=pm,tpar=list(tpar1,tpar2),n.start=350,
                             nrep=1e5,nverb=5000,iseed=c(42,17,69))
    
    # Identical to X1.straushm.trend; tpar given in vector form.
    tpar <- c(2,2,1,tpar1,tpar2)
    X2.straushm.trend <- rmh("straushm",par=par6,w=w,ntypes=2,
                             ptypes=pm,tpar=tpar,n.start=350,
                             nrep=1e5,nverb=5000,iseed=c(42,17,69))
    par7 <- c(3600,0.08)
    w    <- c(0,1,0,1)
    X.dig1 <- rmh("dig1",par=par7,w=w,n.start=300,nrep=1e5,nverb=5000)
    
    par8 <- c(1800,3,0.02,0.04)
    X.dig2 <- rmh("dig2",par=par8,w=w,n.start=300,nrep=1e5,nverb=5000)
    
    par9 <- c(1.25,1.6,0.2,4.5)
    w    <- c(0,10,0,10)
    X1.geyer <- rmh("geyer",par=par9,w=w,n.start=200,nrep=1e5,nverb=5000)

    # Same as a Strauss process with parameters (2.25,0.16,0.7).
    par10 <- c(2.25,0.4,0.7,10000)
    X2.geyer <- rmh("geyer",par=par10,w=w,n.start=70,nrep=1e5,nverb=5000)

    # Starting from the redwood data set, and simulating on a torus:
    par11 <- c(8.1,2.2,0.08,3)
    data(redwood)
    X3.geyer <- rmh("geyer",par=par11,x.start=redwood,periodic=TRUE,
                    nrep=1e5,nverb=5000)

    # As above, conditioning on n:
    X4.geyer <- rmh("geyer",par=par11,x.start=redwood,periodic=TRUE,
                    p=1,nrep=1e5,nverb=5000)
   <testonly>par11 <- c(2,0.2,0.7)
    w     <- c(0,10,0,10)
    X1.strauss <- rmh("strauss",par=par11,w=w,n.start=80,
                      nrep=100,nverb=50)

    # Conditioning on n = 80:
    X2.strauss <- rmh("strauss",par=par11,w=w,n.start=80,
                      p=1,nrep=100,nverb=50)
    
    par12 <- c(2000,0.6,0.07)
    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)
    w     <- owin(poly=list(x=x,y=y))
    X3.strauss <- rmh("strauss",par=par12,w=w,n.start=90,
                      nrep=100,nverb=50)
    
    # Pure hardcore:
    par13 <- c(2,0,0.7)
    w     <- c(0,10,0,10)
    X4.strauss <- rmh("strauss",par=par13,w=w,n.start=60,
                      nrep=100,nverb=50,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))
    X5.strauss <- rmh("strauss",par=par12,w=w,x.start=xxx,
                      nrep=100,nverb=50)
    
    par21 <- c(2,0.2,0.7,0.3)
    w     <- c(0,10,0,10)
    X1.straush <- rmh("straush",par=par21,w=w,n.start=70,
                      nrep=100,nverb=50)
    
    par22 <- c(80,0.36,45,2.5)
    w     <- c(0,250,0,250)
    X2.straush <- rmh("straush",par=par22,w=w,n.start=160,
                      nrep=100,nverb=50)
    
    # Pure hardcore (identical to X3.strauss).
    par23 <- c(2,1,1,0.7)
    w     <- c(0,10,0,10)
    X3.straush <- rmh("straush",par=par23,w=w,n.start=60,
                      nrep=100,nverb=50,iseed=c(42,17,69))
    
    par3 <- c(0.8,0.1,0.5)
    w    <- c(0,10,0,10)
    X.sftcr <- rmh("sftcr",par=par3,w=w,n.start=70,nrep=100,
                   nverb=50)
    
    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)
    par4 <- list(beta=beta,gamma=gmma,r=r)
    w    <- c(0,250,0,250)
    pm   <- c(0.75,0.25)
    X1.straussm <- rmh("straussm",par=par4,w=w,ntypes=2,
                      ptypes=pm,n.start=80,nrep=100,nverb=50)
    
    # Conditioning upon the total number of points being 80:
    X2.straussm <- rmh("straussm",par=par4,w=w,ntypes=2,
                      ptypes=pm,n.start=80,p=1,nrep=100,nverb=50)
    
    # Conditioning upon the number of points of type 1 being 60
    # and the number of points of type 2 being 20:
    X3.straussm <- rmh("straussm",par=par4,w=w,ntypes=2,n.start=c(60,20),
                       fixall=TRUE,nrep=100,p=1,nverb=50)
    
    rhc  <- matrix(c(9.1,5.0,5.0,2.5),2,2)
    par5 <- list(beta=beta,gamma=gmma,r=r,rhc=rhc)
    X.straushm <- rmh("straushm",par=par5,w=w,ntypes=2,
                      ptypes=pm,n.start=80,nrep=100,nverb=50)
    
    beta  <- c(0.0027,0.08)
    par6  <- list(beta=beta,gamma=gmma,r=r,rhc=rhc)
    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.
    w     <- c(0,250,0,250)
    pm    <- c(0.75,0.25)
    X1.straushm.trend <- rmh("straushm",par=par6,w=w,ntypes=2,
                             ptypes=pm,tpar=list(tpar1,tpar2),n.start=350,
                             nrep=100,nverb=50,iseed=c(42,17,69))
    
    # Identical to X1.straushm.trend; tpar given in vector form.
    tpar <- c(2,2,1,tpar1,tpar2)
    X2.straushm.trend <- rmh("straushm",par=par6,w=w,ntypes=2,
                             ptypes=pm,tpar=tpar,n.start=350,
                             nrep=100,nverb=50,iseed=c(42,17,69))
    par7 <- c(3600,0.08)
    w    <- c(0,1,0,1)
    X.dig1 <- rmh("dig1",par=par7,w=w,n.start=300,nrep=100,nverb=50)
    
    par8 <- c(1800,3,0.02,0.04)
    X.dig2 <- rmh("dig2",par=par8,w=w,n.start=300,nrep=100,nverb=50)
    
    par9 <- c(1.25,1.6,0.2,4.5)
    w    <- c(0,10,0,10)
    X1.geyer <- rmh("geyer",par=par9,w=w,n.start=200,nrep=100,nverb=50)

    # Same as a Strauss process with parameters (2.25,0.16,0.7).
    par10 <- c(2.25,0.4,0.7,10000)
    X2.geyer <- rmh("geyer",par=par10,w=w,n.start=70,nrep=100,nverb=50)

    # Starting from the redwood data set, and simulating on a torus:
    par11 <- c(8.1,2.2,0.08,3)
    data(redwood)
    X3.geyer <- rmh("geyer",par=par11,x.start=redwood,periodic=TRUE,
                    nrep=100,nverb=50)

    # As above, conditioning on n:
    X4.geyer <- rmh("geyer",par=par11,x.start=redwood,periodic=TRUE,
                    p=1,nrep=100,nverb=50)</testonly>

Run the code above in your browser using DataLab