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,ntypes=0,ptypes=NULL,tpar=NULL,n.start,expand=NULL,
        periodic=FALSE,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. This must be in a form which can be coerced to an object of class owin by as.owin(). Note that for non-rectangu

ntypes

The number of distinct types to be used if the simulated pattern is to be a multitype 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 shou

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 ter

n.start

The number of ``initial'' points to be randomly (uniformly) generated in the window w. This set of uniformly generated points gives the Metropolis-Hastings algorithm an initial state from which to start. (Actually the number

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 w i

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.

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

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, nrep, p, q, expand, periodic, and iseed.

Warnings

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

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)
    
    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))
    X2.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)
    X3.strauss <- rmh("strauss",par=par13,w=w,n.start=60,
                      nrep=1e5,nverb=5000,iseed=c(42,17,69))
    
    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)
    X.straussm <- rmh("straussm",par=par4,w=w,ntypes=2,
                      ptypes=pm,n.start=80,nrep=1e5,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)
   <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)
    
    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))
    X2.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)
    X3.strauss <- rmh("strauss",par=par13,w=w,n.start=60,
                      nrep=100,nverb=50,iseed=c(42,17,69))
    
    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)
    X.straussm <- rmh("straussm",par=par4,w=w,ntypes=2,
                      ptypes=pm,n.start=80,nrep=100,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)</testonly>

Run the code above in your browser using DataCamp Workspace