rmh.default {spatstat} | R Documentation |
Generates a random point pattern, simulated from a chosen point process model, using the Metropolis-Hastings algorithm.
## Default S3 method: rmh(model, start, control, verbose=TRUE, track=FALSE, ...)
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 |
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:
A character string specifying the choice of interpoint interaction for the point process.
Parameter values for the conditional intensity function.
(Optional) window in which the pattern is
to be generated. An object of class "owin"
,
or data acceptable to as.owin
.
Data specifying the spatial trend in the model, if it has a trend.
This may be a function, a pixel image (of class "im"
),
(or a list of functions or images if the model
is multitype).
If the trend is a function or functions,
any auxiliary arguments ...
to rmh.default
will be passed to these functions, which
should be of the form function(x, y, ...)
.
List of possible types, for a multitype point process.
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:
Number of points in the initial point pattern.
A single integer, or a vector of integers giving the
numbers of points of each type in a multitype point pattern.
Incompatible with x.start
.
Initial point pattern configuration.
Incompatible with n.start
.
x.start
may be a point pattern (an
object of class "ppp"
), or data which can be coerced
to this class by as.ppp
, or an object with
components x
and y
, or a two-column matrix.
In the last two cases, the window for the pattern is determined
by model$w
.
In the first two cases, if model$w
is also present,
then the final simulated pattern will be clipped to
the window model$w
.
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:
The probability of proposing a “shift” (as opposed to a birth or death) in the Metropolis-Hastings algorithm.
The conditional probability of proposing a death (rather than a birth) given that birth/death has been chosen over shift.
The number of repetitions or iterations to be made by the Metropolis-Hastings algorithm. It should be large.
Either a numerical expansion factor, or
a window (object of class "owin"
). Indicates that
the process is to be simulated on a larger domain than the
original data window w
, then clipped to w
when the algorithm has finished.
The default is to expand the simulation window if the model is stationary and non-Poisson (i.e. it has no trend and the interaction is not Poisson) and not to expand in all other cases.
If the model has a trend, then in order for expansion to be feasible, the trend must be given either as a function, or an image whose bounding box is large enough to contain the expanded window.
A logical scalar; if periodic
is TRUE
we simulate a process on the torus formed by identifying
opposite edges of a rectangular window.
A vector of probabilities (summing to 1) to be used in assigning a random type to a new point.
A logical scalar specifying whether to condition on the number of points of each type.
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.
If this argument is present, then
conditional simulation will be performed, and x.cond
specifies the conditioning points and the type of conditioning.
For full details of these parameters, see rmhcontrol
.
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.
There are several kinds of conditional simulation.
Simulation conditional upon the number of points,
that is, holding the number of points fixed.
To do this, set control$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 setting start$n.start
to be a
scalar, or by setting the initial pattern start$x.start
.
In the case of multitype processes, it is possible to simulate the
model conditionally upon the number of points of each type,
i.e. holding the number of points of each type
to be fixed. To do this, set control$p
equal to 1
and control$fixall
to be TRUE
.
The number of points is then determined by the starting state, which
may be specified either by setting start$n.start
to be an
integer vector, or by setting the initial pattern start$x.start
.
Simulation
conditional 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, set control$x.cond
to equal the
specified point pattern y, making sure that it is an object of class
"ppp"
and that the window as.owin(control$x.cond)
is the conditioning window V.
Simulation conditional 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, set control$x.cond
to be a
data.frame
containing 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.
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.
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.
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
.
Adrian Baddeley Adrian.Baddeley@csiro.au http://www.maths.uwa.edu.au/~adrian/ and Rolf Turner r.turner@auckland.ac.nz
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.
rmh
,
rmh.ppm
,
rStrauss
,
ppp
,
ppm
,
AreaInter
,
BadGey
,
DiggleGatesStibbard
,
DiggleGratton
,
Fiksel
,
Geyer
,
Hardcore
,
LennardJones
,
MultiHard
,
MultiStrauss
,
MultiStraussHard
,
PairPiece
,
Poisson
,
Softcore
,
Strauss
,
StraussHard
,
Triplets
## Not run: nr <- 1e5 nv <- 5000 ns <- 200 ## End(Not run) 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)) # 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)) # Triplets process modtrip <- list(cif="triplets",par=list(beta=2,gamma=0.2,r=0.7), w=c(0,10,0,10)) X.triplets <- rmh(model=modtrip, 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. ## Not run: 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. ## End(Not run)