rmh.default {spatstat} | R Documentation |
Generates a random point pattern, simulated from a chosen point process model, using the Metropolis-Hastings algorithm.
rmh.default(model,start,control,...)
model |
List of parameters specifying the point process model
that is to be simulated:
|
start |
List of parameters determining the initial state of
the algorithm:
n.start and x.start are incompatible.
See Details for details.
|
control |
List of parameters controlling the iterative behaviour
and termination of the algorithm:
|
... |
Further arguments; currently ignored. |
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:
'strauss'
'straush'
'sftcr'
'straussm'
'straushm'
'dgs'
'diggra'
'geyer'
See the section Extensions for the possibility of extending this list of options.
par
should be
as follows, for each of the available conditional intensity
functions:
beta,gamma,r
which are respectively the ``base''
intensity, the pair-wise interaction parameter and the
interaction radius. Note that gamma
must be less than
or equal to 1.
beta,gamma,r,hc
where beta
, gamma
,
and r
are as for the Strauss process, and hc
is
the hardcore radius. Of course hc
must be less than
r
.
beta,sigma,kappa
. Again
beta
is a ``base'' intensity. The pairwise
interaction between two points u != v is
-(sigma/||u-v||)^(2/kappa)
Note that it is necessary that 0 < kappa <1.
beta
:
A vector of ``base'' intensities, one for each possible
type.
gamma
:
A symmetric matrix of interaction parameters,
with gamma_ij pertaining to the interaction between
type i and type j.
radii
:
A symmetric matrix of interaction radii, with
entries r_ij pertaining to the interaction between type
i and type j.
beta
and gamma
as for straussm
and
two ``radii'' components:
iradii
: the interaction radii
hradii
: the hardcore radii
which are both symmetric matrices of nonnegative numbers.
The entries of hradii
must be less than the
corresponding entries
of iradii
.
beta
and rho
. This process has pairwise interaction
function
e(t) = sin^2((pi * t)/(2 * rho))
for t < rho, and equal to 1 for t >= rho.
beta
,
kappa
, delta
and rho
. This process has
pairwise interaction function e(t) equal to 0
for t < delta, equal to
((t-delta)/(rho-delta))^kappa
for delta <= t < rho, and equal to 1 for delta >= rho. Note that here we use the symbol kappa where Diggle, Gates, and Stibbard use beta since we reserve the symbol beta for an intensity parameter.
beta
, gamma
, r
, and sat
.
The components beta
, gamma
, r
are as for
the Strauss model, and sat
is the ``saturation''
parameter. The model is Geyer's ``saturation'' point process
model. It is ``like a Strauss model, but with an upper bound
to the number of r-close neighbors of any point.''
Explicitly, a saturation point process with interaction radius r, saturation threshold s, and parameters beta and gamma, is the point process in which each point x[i] in the pattern X contributes a factor
beta gamma^min(s,t(x[i],X))
to the probability density of the point pattern, where t(x[i],X) denotes the number of ``r-close neighbours'' of x[i] in the pattern X.
tpar
. If the user wishes to omit some of the terms
in such a polynomial, the corresponding coefficients should be
specified as 0. For example, to specify a trend 2x the
trend parmeters should be specified as c(2,0)
.
If the process is unmarked tpar
must be a vector
consisting solely of the coefficients of the trend polynomial in
the correct order.
If the process is multitype, so that tpar
is given as
a list, the i-th component of this list must
consist of the coefficients for the trend corresponding to the
i-th type. If there is to be no trend for the
i-th type, then the i-th component
of the list must be NULL (but must be present in the list as
the NULL object.)
owin
by as.owin()
. Note that for non-rectangular
windows we cannot (as yet) generate a pattern from a process
which exists only within that window. The domain of the
theoretical process must (at present) be at least as large as
the ``enclosing box'' of the window.
The argument start
determines the initial state of the
Metropolis-Hastings algorithm. Possible components are
w
.
For a multitype point process, n.start
may be a vector (of
length equal to the number of types)
giving the number of points of each type to be generated.
Incompatible with x.start
.
n.start
is the number of ``initial'' points to be randomly
(uniformly) generated in the window w
, or a vector (of
length equal to the number of types for a multitype process)
giving the number of points of each type to be generated. These
uniformly generated points form an initial state or configuration
for the Metropolis-Hastings algorithm. Precisely one of
n.start
or x.start
must be specified.
A vector-valued n.start
is meaningful only
if p
(the probability of a shift as opposed to a birth or
death) is equal to 1 (so that we are conditioning on the number
of points). When p < 1, if n.start
is vector valued
then it is effectively replaced by its sum.
The resulting set of uniformly generated points gives the
Metropolis-Hastings algorithm an initial state from which to
start. (Actually, when p < 1, the number n.start
gets multiplied by the ratio of the area of the enclosing box for
the window to the area of the window, and then by the factor
expand
. Then that many points are uniformly generated in
the expanded window; see below.) The value of
n.start
should be roughly equal to (an educated guess at)
the expected number of points which will be generated inside the
window.
"ppp"
,
or data which can be coerced to this class by
as.ppp
).
Incompatible with n.start
.
x.start
is a point pattern (an object of class ppp
,
or an object which can be coerced to this class by as.ppp()
).
This object provides an alternative means of specifying the initial
``state'' or configuration for the Metropolis-Hastings algorithm.
If x.start
is specified, but x.start$window
is NULL
,
then this gap is filled in by the component w
of start
.
If x.start$window
is present, and if w
is specified
as well, the latter is used as a window to which to clip the final
simulated pattern. Thus in such circumstances it is only
sensible to specify a value of w
which is contained in
x.start$window
. However no checking is done for this.
The simulated pattern is constructed in the enclosing rectangle
of x.start$window
. No expansion takes place. (The
argument expand
is forced to equal 1.) As indicated
above, at the end of the simulation, the resulting pattern is
clipped to the window w
if this is given (or to
x.start$window
if this differs from its enclosing
rectangle).
sample
.
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:
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 1, then we
are simulating the latter (unless periodic
is TRUE
;
see below). The larger expand
is, the better we
approximate the former. Note that any value of expand
smaller than 1 is treated as if it were 1.
The area of the expanded window is equal to expand
times
the area of the enclosing box; width and height are stretched
proportionately. Points are generated by the Metropolis-Hastings
algorithm in the expanded window, and then ``clipped'' down to
the original window w
(or to x.start$window
) when the
algorithm has finished. The argument expand
defaults to 2
if periodic
is FALSE
and p < 1 and to 1 if
periodic
is TRUE
or if p < 1, or if the
starting configuration is specified via x.start
. Trying to
set expand
greater than 1 when periodic
is
TRUE
or p = 1 generates an error. A specified value
of expand
is simply ignored if x.start
is given.
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 (i.e. x.start$window
if
x.start
is specified; otherwise w
) is not
rectangular, an error is given.
fixall=TRUE
. In this case, n.start
must be a
vector whose entries are these numbers.
We can only condition on the number of points if the simulation
takes place in the original window (as opposed to taking place in
a larger window and then being clipped to the original). Hence,
if p = 1 then expand
defaults to 1 and it is an
error to specify a value of expand
which is greater than
1. Likewise, if p = 1, it is (currently) an error to specify
anything other than a rectangular window.
p
is equal to 1.
ptypes
is close to the
relative frequencies of the types which will result from the
simulation.
fixall
is set equal to TRUE
when it is
not meaningful.
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
.
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
expand
must be equal to 1; it
defaults to 1, and it is an error to specify a value larger
than 1.}
n.start
.
}
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.
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
.)
Adrian Baddeley adrian@maths.uwa.edu.au http://www.maths.uwa.edu.au/~adrian/ and Rolf Turner rolf@math.unb.ca http://www.math.unb.ca/~rolf
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{o}ller, J. (1994) Simulation procedures and likelihood inference for spatial point processes. Scandinavian Journal of Statistics 21, 359373.
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 79140.
rmh
,
rmh.ppm
,
ppp
,
mpl
,
Strauss
,
Softcore
,
StraussHard
,
MultiStrauss
,
MultiStraussHard
,
DiggleGratton
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))