rmh.ppm
Simulate from a Fitted Point Process Model
Given a point process model fitted to data, generate a random simulation of the model, using the Metropolis-Hastings algorithm.
- Keywords
- spatial
Usage
rmh.ppm(model,start=NULL,control=NULL,..., verbose=TRUE, project=TRUE)Arguments
- model
- A fitted point process model (object of class
"ppm", seeppm.object) which it is desired to simulate. This fitted model is usually the result of a call to - start
- A list of arguments determining the initial state
of the Metropolis-Hastings algorithm. See
rmh.defaultfor description of these arguments. Defaults to list(x.start=model$Q$data) - control
- A list of arguments controlling the running of
the Metropolis-Hastings algorithm. See
rmh.defaultfor description of these arguments. - ...
- Further arguments passed to
rmh.default. - verbose
- Logical flag indicating whether to print progress reports.
- project
- Logical flag indicating what to do if the fitted model is
invalid (in the sense that the values of the fitted coefficients do not
specify a valid point process).
If
project=TRUEthe closest valid model will be simulated; if
Details
This function generates simulated realisations from a point
process model that has been fitted to point pattern data. It is
a method for the generic function rmh for the
class "ppm" of fitted point process models. To simulate
other kinds of point process models, see rmh
or rmh.default.
The argument model describes the fitted model. It must be
an object of class "ppm" (see ppm.object),
and will typically be the result of a call to the point process
model fitting function mpl.
The current implementation enables simulation from any fitted model
involving the interactions
DiggleGratton,
Geyer,
MultiStrauss,
MultiStraussHard,
Poisson,
Strauss,
StraussHard
and Softcore,
including nonstationary models, provided the trend does not
involve spatial covariates. See the examples.
It is possible that the fitted coefficients of a point process model
may be ``illegal'', i.e. that there may not exist a
mathematically well-defined point process with the given parameter
values. For example, a Strauss process with interaction
parameter $\gamma > 1$ does not exist,
but the model-fitting procedure used in mpl will sometimes
produce values of $\gamma$ greater than 1.
In such cases, if project=FALSE then an error will occur,
while if project=TRUE then rmh.ppm will find
the nearest legal model and simulate
this model instead. (The nearest legal model is obtained by
projecting the vector of coefficients onto the set of
valid coefficient vectors. The result is usually the Poisson process
with the same fitted intensity.)
The arguments start and control are lists of
parameters determining the initial state and the iterative
behaviour, respectively, of the Metropolis-Hastings algorithm.
They are passed directly to rmh.default.
See rmh.default for details of these parameters.
Note that if the model has a trend then the expand
component of control will be ignored.
The value of expand must be equal to 1 in this setting.
After extracting the relevant information from the fitted model
object model, rmh.ppm simply invokes the default
rmh algorithm rmh.default, unless the model
is Poisson. If the model is Poisson then rmh.default is
not needed, and a stand-in function pseudo.rmh is invoked
to do the simple simulation that is required.
See rmh.default for further information about the
implementation, or about the Metropolis-Hastings algorithm.
Value
- A point pattern (an object of class
"ppp"; seeppp.object).
Warnings
See Warnings in rmh.default.
See Also
rmh,
rmh.default,
ppp.object,
mpl,
Poisson,
Strauss,
StraussHard,
Softcore,
Geyer,
DiggleGratton
Examples
require(spatstat)
data(swedishpines)
X <- swedishpines
plot(X, main="Swedish Pines data")
# Poisson process
fit <- mpl(X, ~1, Poisson())
Xsim <- rmh(fit)
plot(Xsim, main="simulation from fitted Poisson model")
# Strauss process
fit <- mpl(X, ~1, Strauss(r=7), rbord=7)
Xsim <- rmh(fit, start=list(n.start=X$n), control=list(nrep=1e3))
plot(Xsim, main="simulation from fitted Strauss model")
# Strauss - hard core process
fit <- mpl(X, ~1, StraussHard(r=7,hc=2), rbord=7)
Xsim <- rmh(fit, start=list(n.start=X$n), control=list(nrep=1e3))
plot(Xsim, main="simulation from fitted Strauss hard core model")
# Geyer saturation process
fit <- mpl(X, ~1, Geyer(r=7,sat=2), rbord=7)
Xsim <- rmh(fit, start=list(n.start=X$n), control=list(nrep=1e3))
plot(Xsim, main="simulation from fitted Geyer model")
# soft core interaction process
fit <- mpl(X, ~1, Softcore(kappa=0.1))
Xsim <- rmh(fit, start=list(n.start=X$n), control=list(nrep=1e3))
plot(Xsim, main="simulation from fitted Soft Core model")
data(cells)
plot(cells)
# Diggle-Gratton pairwise interaction model
fit <- mpl(cells, ~1, DiggleGratton(0.05, 0.1))
Xsim <- rmh(fit, start=list(n.start=cells$n), control=list(nrep=1e3))
plot(Xsim, main="simulation from fitted Diggle-Gratton model")
# marked point pattern
data(amacrine)
Y <- amacrine
plot(Y, main="Amacrine data")
# marked Poisson models
fit <- mpl(Y)
Ysim <- rmh(fit)
plot(Ysim, main="simulation from mpl(Y)")
fit <- mpl(Y,~marks)
Ysim <- rmh(fit)
plot(Ysim, main="simulation from mpl(Y, ~marks)")
fit <- mpl(Y,~polynom(x,y,2))
Ysim <- rmh(fit)
plot(Ysim, main="simulation from mpl(Y, ~polynom(x,y,2))")
fit <- mpl(Y,~marks+polynom(x,y,2))
Ysim <- rmh(fit)
plot(Ysim, main="simulation from mpl(Y, ~marks+polynom(x,y,2))")
fit <- mpl(Y,~marks*polynom(x,y,2))
Ysim <- rmh(fit)
plot(Ysim, main="simulation from mpl(Y, ~marks*polynom(x,y,2))")
# multitype Strauss models
MS <- MultiStrauss(types = levels(Y$marks),
radii=matrix(0.07, ncol=2, nrow=2))
fit <- mpl(Y, ~marks, MS)
Ysim <- rmh(fit, start=list(n.start=X$n), control=list(nrep=1e3))
plot(Ysim, main="simulation from fitted Multitype Strauss")
fit <- mpl(Y,~marks*polynom(x,y,2), MS)
Ysim <- rmh(fit, start=list(n.start=X$n), control=list(nrep=1e3))
plot(Ysim, main="simulation from fitted inhomogeneous Multitype Strauss")