| rmh {spatstat} | R Documentation |
Simulates point patterns from a (still small) range of point process models.
rmh(cif,par,w,ntypes=0,ptypes=NULL,tpar=NULL,n.start,expand=NULL,
periodic=F,nrep=1e6,p=0.9,q=0.5,iseed=NULL,nverb=0)
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),
straush (Strauss process with hardcore), sftcr
(Softcore process), straussm (multitype Strauss process),
straushm (multitype Strauss process with hardcore),
dig1, dig2 (see the 2nd of the References),
and geyer (see the 3rd of the References).
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 or ratfor code for the new conditional intensity function, the user would have to modify the code in the files cif.r and rmh.R
appropriately. (And recreate/recompile the dynamically loadable
shared object spatstat.so.)
|
par |
A set of parameter values appropriate to the conditional intensity function being invoked. (See Details.) |
w |
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-rectangular windows we cannot (as yet)
generated 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.
|
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 should be improved if ptypes
is close to the relative frequencies of the types which will
result from the simulation.
|
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
terms x, y, x^2, xy, y^2, x^3, x^2y, xy^2, y^3, ..., y^n
where n is the degree of the trend. (Thus the
possible lengths of a sequence of coefficients
are 0, 2, 5, 9, 14, ...).
If the process is unmarked tpar must be a vector
consisting solely of the coefficients of the trend polynomial
in the correct order. The degree of the polynomial will be
calculated from the length of tpar.
If the process is multitype and if tpar is given as a vector,
then ntypes and and the degrees of each of the
polynomial should comprise the first 1+ntypes entries
of tpar. After that, the coefficients corresponding to
each mark should be catenated successively. Note that some
of the degrees are allowed to be 0.
If tpar is given as a list, the i-th
component of this list should 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 should be NULL (but must
be present in the list as the NULL object.)
|
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 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.
|
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 1,
then we are simulating the latter (unless periodic
in 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 when the algorithm has finished.
The argument expand defaults to 2 when periodic
is FALSE and to 1 when periodic is TRUE. (Trying
to set expand greater than 1 when periodic is TRUE
generates an error.)
|
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 is not rectangular, an error is given.
|
nrep |
The number of repetitions or steps (changes of state) to be made by the Metropolis-Hastings algorithm. It should be large. |
p |
The probability of of proposing a ``shift'' (as opposed to a birth or death) in the Metropolis-Hastings algorithm. |
q |
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. |
The parameter vector (or list) par should be
as follows, for each of the available conditional intensity
functions:
-(sigma/||u-v||)^(2/kappa)
Note that it is necessary that 0 < kappa <1.
par is
best given as a list with components
If par is to be given as a vector, then this vector should
consist of the vector beta, catenated with the upper
triangle of the matrix gamma strung out in row
order, and finally catenated with the upper triangle of the
matrix r likewise strung out in row order.
par are much as for straussm except
that there is an extra component rhc which is the matrix
of hardcore radii. Note that in the vector form, the hardcore
radii must be catenated after the interaction radii.
e(t) = sin^2((pi t)/(2 rho))
for t < rho, and equal to 1 for t >= rho. The parameters are beta and rho.
((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.
The complete parameter set is beta, kappa, delta and rho.
More 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.
A list of class ppp with the usual components
window |
The window w in the form of
an object of class owin.
|
n |
The number of generated points (in the window w). |
x |
The x-coordinates of the generated points. |
y |
The y-coordinates of the generated points. |
marks |
The 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.
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.
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, 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
ppp,
mpl,
Strauss,
Softcore,
StraussHard,
MultiStrauss,
MultiStraussHard
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)