Sets up a list of parameters controlling the iterative behaviour of the Metropolis-Hastings algorithm.
rmhcontrol(…) # S3 method for default
rmhcontrol(…, p=0.9, q=0.5, nrep=5e5,
expand=NULL, periodic=NULL, ptypes=NULL,
x.cond=NULL, fixall=FALSE, nverb=0,
nsave=NULL, nburn=nsave, track=FALSE,
pstage=c("block", "start"))
Arguments passed to methods.
Probability of proposing a shift (as against a birth/death).
Conditional probability of proposing a death given that a birth or death will be proposed.
Total number of steps (proposals) of Metropolis-Hastings algorithm that should be run.
Simulation window or expansion rule.
Either a window (object of class "owin")
or a numerical expansion factor, specifying that
simulations are to be performed in a domain other than the
original data window, then clipped to the original data window.
This argument is passed to rmhexpand.
A numerical expansion factor can be in several formats:
see rmhexpand.
Logical value (or NULL) indicating whether to simulate
``periodically'', i.e. identifying opposite edges of the rectangular
simulation window. A NULL value means ``undecided.''
For multitype point processes, the distribution of the mark attached to a new random point (when a birth is proposed)
Conditioning points for conditional simulation.
(Logical) for multitype point processes, whether to fix the number of points of each type.
Progress reports will be printed every nverb
iterations
If these values are specified, then
intermediate states of the simulation algorithm will be saved
every nsave iterations, after an initial burn-in period of
nburn iterations.
Logical flag indicating whether to save the transition history of the simulations.
Character string specifying when to generate
proposal points. Either "start" or "block".
An object of class "rmhcontrol", which is essentially
a list of parameter values for the algorithm.
There is a print method for this class, which prints
a sensible description of the parameters chosen.
For a Gibbs point process \(X\), the Metropolis-Hastings algorithm easily accommodates several kinds of conditional simulation:
We fix the total number of points \(N(X)\) to be equal to \(n\). We simulate from the conditional distribution of \(X\) given \(N(X) = n\).
In a multitype point process, where \(Y_j\) denotes the process of points of type \(j\), we fix the number \(N(Y_j)\) of points of type \(j\) to be equal to \(n_j\), for \(j=1,2,\ldots,m\). We simulate from the conditional distribution of \(X\) given \(N(Y_j)=n_j\) for \(j=1,2,\ldots,m\).
We require that the point process \(X\) should, within a specified sub-window \(V\), coincide with a specified point pattern \(y\). We simulate from the conditional distribution of \(X\) given \(X \cap V = y\).
We require that the point process \(X\) include a specified list of points \(y\). We simulate from the point process with probability density \(g(x) = c f(x \cup y)\) where \(f\) is the probability density of the original process \(X\), and \(c\) is a normalising constant.
To achieve each of these types of conditioning we do as follows:
Set p=1.
The number of points is determined by the initial state
of the simulation: see rmhstart.
Set p=1 and fixall=TRUE.
The number of points of each type is determined by the initial state
of the simulation: see rmhstart.
Set x.cond to be a point pattern (object of
class "ppp"). Its window V=Window(x.cond) becomes the
conditioning subwindow \(V\).
Set x.cond to be a list(x,y) or data.frame
with two columns containing the coordinates of the points, or a
list(x,y,marks) or data.frame with three columns
containing the coordinates and marks of the points.
The arguments x.cond, p and fixall can be
combined.
The Metropolis-Hastings algorithm, implemented as rmh,
generates simulated realisations of point process models.
The function rmhcontrol
sets up a list of parameters which control the
iterative behaviour
and termination of the Metropolis-Hastings algorithm, for use in a
subsequent call to rmh. It also checks that the
parameters are valid.
(A separate function rmhstart
determines the initial state of the algorithm,
and rmhmodel determines the model to be simulated.)
The parameters are as follows:
The probability of proposing a ``shift'' (as opposed to a birth or death) in the Metropolis-Hastings algorithm.
If \(p = 1\) then the algorithm only alters existing points,
so the number of points never changes, i.e. we are
simulating conditionally upon the number of points.
The number of points is determined by the initial state
(specified by rmhstart).
If \(p=1\) and fixall=TRUE and the model
is a multitype point process model, then the algorithm
only shifts the locations of existing points and does not
alter their marks (types).
This is equivalent to simulating conditionally
upon the number of points of each type.
These numbers are again specified by the initial state.
If \(p = 1\) then no expansion of the simulation window
is allowed (see expand below).
The default value of p can be changed by setting
the parameter rmh.p in spatstat.options.
The conditional probability of proposing a death
(rather than a birth) given that a shift is not proposed.
This is of course ignored if p is equal to 1.
The default value of q can be changed by setting
the parameter rmh.q in spatstat.options.
The number of repetitions or iterations to be made by the Metropolis-Hastings algorithm. It should be large.
The default value of nrep can be changed by setting
the parameter rmh.nrep in spatstat.options.
Either a number or a window (object of class "owin").
Indicates that the process is to be simulated on a
domain other than the original data window w,
then clipped to w when the algorithm has finished.
This would often be done in order to approximate the
simulation of a stationary process (Geyer, 1999)
or more generally a process existing in the
whole plane, rather than just in the window w.
If expand is a window object, it is taken as the
larger domain in which simulation is performed.
If expand is numeric, it is interpreted
as an expansion factor or expansion distance
for determining the simulation domain from the data window.
It should be a named scalar, such as
expand=c(area=2), expand=c(distance=0.1),
expand=c(length=1.2). See rmhexpand() for
more details. If the name is omitted, it defaults to area.
Expansion is not permitted if the number of points has been
fixed by setting p = 1 or if the
starting configuration has been specified via the
argument x.start in rmhstart.
If expand is NULL, this is interpreted to mean
“not yet decided”. An expansion rule will be determined
at a later stage, using appropriate defaults.
See rmhexpand.
A logical value (or NULL)
determining whether to simulate “periodically”.
If periodic is TRUE, and if the simulation window
is a rectangle, then the simulation algorithm effectively
identifies opposite edges of the rectangle. Points
near the right-hand edge of the rectangle are deemed to be close
to points near the left-hand edge. Periodic simulation usually
gives a better approximation to a stationary point process.
For periodic simulation, the simulation window must be a rectangle.
(The simulation window is determined by expand as described
above.)
The value NULL means ‘undecided’.
The decision is postponed until rmh is called.
Depending on the point process model to be simulated,
rmh will then set periodic=TRUE if the simulation window
is expanded and the expanded simulation window is rectangular;
otherwise periodic=FALSE.
Note that periodic=TRUE is only permitted when the
simulation window (i.e. the expanded window) is rectangular.
A vector of probabilities (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/nt\) where \(nt\) is the number
of types for the process. Convergence of the simulation
algorithm should be improved if ptypes is close to the
relative frequencies of the types which will result from the
simulation.
If this argument is given,
then conditional simulation will be performed,
and x.cond specifies the location of the
fixed points as well as the type of conditioning.
It should be either a point pattern
(object of class "ppp") or a list(x,y)
or a data.frame.
See the section on Conditional Simulation.
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 when it is
not meaningful.
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 these integers are given, then the
current state of the simulation algorithm (i.e. the current
random point pattern) will be saved every nsave iterations,
starting from iteration nburn.
Logical flag indicating whether to save the transition history of the simulations (i.e. information specifying what type of proposal was made, and whether it was accepted or rejected, for each iteration).
Character string specifying the stage of the algorithm
at which the randomised proposal points should be generated.
If pstage="start" or if nsave=0,
the entire sequence of nrep
random proposal points is generated at the start of the
algorithm. This is the original
behaviour of the code, and should be used in order to maintain
consistency with older versions of spatstat.
If pstage="block" and nsave > 0, then
a set of nsave random proposal points will be generated
before each block of nsave iterations. This is much more
efficient.
The default is pstage="block".
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.
# NOT RUN {
# parameters given as named arguments
c1 <- rmhcontrol(p=0.3,periodic=TRUE,nrep=1e6,nverb=1e5)
# parameters given as a list
liz <- list(p=0.9, nrep=1e4)
c2 <- rmhcontrol(liz)
# parameters given in rmhcontrol object
c3 <- rmhcontrol(c1)
# }
Run the code above in your browser using DataLab