spatstat (version 1.9-4)

rmhmodel: Define Point Process Model for Metropolis-Hastings Simulation.

Description

Builds a description of a point process model for use in simulating the model by the Metropolis-Hastings algorithm.

Usage

rmhmodel(model)
   ## S3 method for class 'default':
rmhmodel(\dots, cif=NULL, par=NULL, w=NULL, trend=NULL, types=NULL)

Arguments

model
An existing description of the model in some other format. Incompatible with the arguments listed below.
...
There should be no other arguments.
cif
Character string specifying the choice of model
par
Parameters of the model
w
Spatial window in which to simulate
trend
Specification of the trend in the model
types
A vector of factor levels defining the possible marks, for a multitype process.

Value

  • An object of class "rmhmodel", which is essentially a list of parameter values for the model. There is a print method for this class, which prints a sensible description of the model chosen.

synopsis

rmhmodel(model, ...) ## S3 method for class 'default': rmhmodel(model, \dots, cif=NULL, par=NULL, w=NULL, trend=NULL, types=NULL)

Warnings in Respect of ``lookup''

The syntax of rmh.default in respect of the lookup cif has changed from the previous release of spatstat (versions 1.4-7 to 1.5-1). Read the Details carefully. In particular it is now required that the first entry of the r component of par be strictly positive. (This is the opposite of what was required in the previous release, which was that this first entry had to be 0.)

It is also now required that the entries of r be sorted into ascending order. (In the previous release it was assumed that the entries of r and h were in corresponding order and the two vectors were sorted commensurately. It was decided that this is dangerous sand unnecessary.)

Note that if you specify the lookup pairwise interaction function via stepfun() the arguments x and y which are passed to stepfun() are slightly different from r and h: length(y) is equal to 1+length(x); the final entry of y must be equal to 1 --- i.e. this value is explicitly supplied by the user rather than getting tacked on internally.

The step function returned by stepfun() must be right continuous (this is the default behaviour of stepfun()) otherwise an error is given.

Details

Simulated realisations of many point process models can be generated using the Metropolis-Hastings algorithm rmh. The algorithm requires the model to be specified in a particular format.

This function rmhmodel creates a description of the point process model in the form required by rmh, and checks that the parameters of the model are valid.

The point process model should be specified either by the first argument model or by the other arguments cif, par etc. If model is a fitted point process model (object of class "ppm" obtained by a call to the model-fitting function ppm) then a description of the point process model will be extracted from this object.

If model is a list, then it should have components named cif, par and optionally w, trend, types with the same interpretation as described below.

The argument cif is a character string specifying the choice of interpoint interaction for the point process. The current options are [object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object]

The argument par supplies parameter values appropriate to the conditional intensity function being invoked. These are: [object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object]

The optional argument trend determines the spatial trend in the model, if it has one. It should be a function or image (or a list of such, if the model is multitype) to provide the value of the trend at an arbitrary point. [object Object],[object Object] Note that the trend or trends must be non-negative; no checking is done for this. The optional argument w specifies the 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.

The optional argument types specifies the possible types in a multitype point process. If the model being simulated is multitype, and types is not specified, then this vector defaults to 1:ntypes where ntypes is the number of types.

References

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. 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, rmhcontrol, rmhstart, ppm, Strauss, Softcore, StraussHard, MultiStrauss, MultiStraussHard, DiggleGratton, PairPiece

Examples

Run this code
# Strauss process:
   mod01 <- rmhmodel(cif="strauss",par=c(beta=2,gamma=0.2,r=0.7),
                 w=c(0,10,0,10))
   # Equivalent to:
   a <- list(cif="strauss",par=c(beta=2,gamma=0.2,r=0.7),
                 w=c(0,10,0,10))
   mod01 <- rmhmodel(a)
   # The above could also be simulated using 'rStrauss'

   # Strauss with hardcore:
   mod04 <- list(cif="straush",par=c(beta=2,gamma=0.2,r=0.7,hc=0.3),
                w=owin(c(0,10),c(0,5)))
   mod04 <- rmhmodel(mod04)

   # Soft core:
   par3 <- c(0.8,0.1,0.5)
   w    <- square(10)
   mod07 <- rmhmodel(cif="sftcr",
                     par=c(beta=0.8,sigma=0.1,kappa=0.5),
                     w=w)
   
   # 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 <- rmhmodel(cif="straussm",
                     par=list(beta=beta,gamma=gmma,radii=r),
                     w=square(250))
   # specify types
   mod09 <- rmhmodel(cif="straussm",
                     par=list(beta=beta,gamma=gmma,radii=r),
                     w=square(250),
                     types=c("A", "B"))

   
   # Multitype Strauss hardcore with trends for each type:
   beta  <- c(0.27,0.08)
   ri    <- matrix(c(45,45,45,45),2,2)
   rhc  <- matrix(c(9.1,5.0,5.0,2.5),2,2)
   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 <- rmhmodel(cif="straushm",par=list(beta=beta,gamma=gmma,
                 iradii=ri,hradii=rhc),w=c(0,250,0,250),
                 trend=list(tr3,tr4))

   # 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 <- rmhmodel(cif="lookup",par=list(beta=4000,h=h,r=r),w=c(0,1,0,1))

Run the code above in your browser using DataCamp Workspace