rmh.ppm {spatstat} | R Documentation |
Given a point process model fitted to data, generate a random simulation of the model, using the Metropolis-Hastings algorithm.
## S3 method for class 'ppm': rmh(model,start=NULL, control=rmhcontrol(), ..., verbose=TRUE, project=TRUE)
model |
A fitted point process model (object of class
"ppm" , see ppm.object ) which it is desired
to simulate. This fitted model is usually the result of a call
to ppm . See Details below.
|
start |
Data determining the initial state
of the Metropolis-Hastings algorithm. See
rmhstart for description of these arguments.
Defaults to list(x.start=data.ppm(model))
|
control |
Data controlling the running of
the Metropolis-Hastings algorithm. See rmhcontrol
for 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=TRUE the closest valid model will be simulated;
if project=FALSE an error will occur.
|
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 ppm
.
The current implementation enables simulation from any fitted model
involving the interactions
DiggleGratton
,
Geyer
,
MultiStrauss
,
MultiStraussHard
,
PairPiece
,
Poisson
,
Strauss
,
StraussHard
and Softcore
,
including nonstationary models. 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 ppm
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 rmhstart
and
rmhcontrol
respectively.
See rmhstart
and
rmhcontrol
for details of these parameters.
Note that if you specify control$expand > 1
(so that the
model will be simulated on a window larger than the original data
window) then the model must be capable of extrapolation to this
larger window. This excludes models which depend on external covariates.
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 the Metropolis-Hastings
algorithm is not needed, and the model is simulated directly, using
one of rpoispp
, rmpoispp
,
rpoint
or rmpoint
.
See rmh.default
for further information about the
implementation, or about the Metropolis-Hastings algorithm.
A point pattern (an object of class "ppp"
; see
ppp.object
).
See Warnings in rmh.default
.
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
rmh
,
rmhmodel
,
rmhcontrol
,
rmhstart
,
rmh.default
,
ppp.object
,
ppm
,
PairPiece
,
Poisson
,
Strauss
,
StraussHard
,
Softcore
,
Geyer
,
DiggleGratton
data(swedishpines) X <- swedishpines plot(X, main="Swedish Pines data") # Poisson process fit <- ppm(X, ~1, Poisson()) Xsim <- rmh(fit) plot(Xsim, main="simulation from fitted Poisson model") # Strauss process fit <- ppm(X, ~1, Strauss(r=7), rbord=7) Xsim <- rmh(fit, control=list(nrep=1e3)) plot(Xsim, main="simulation from fitted Strauss model") ## Not run: # Strauss process simulated on a larger window # then clipped to original window Xsim <- rmh(fit, control=list(nrep=1e3, expand=2, periodic=TRUE)) # Strauss - hard core process fit <- ppm(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 <- ppm(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 Q <- quadscheme(X, nd=50) fit <- ppm(Q, ~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 <- ppm(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") X <- rSSI(0.05, 100) plot(X, main="new data") # piecewise-constant pairwise interaction function fit <- ppm(X, ~1, PairPiece(seq(0.02, 0.1, by=0.01))) Xsim <- rmh(fit, control=list(nrep=1e3)) plot(Xsim, main="simulation from fitted pairwise model") # marked point pattern data(amacrine) Y <- amacrine plot(Y, main="Amacrine data") # marked Poisson models fit <- ppm(Y) Ysim <- rmh(fit) plot(Ysim, main="simulation from ppm(Y)") fit <- ppm(Y,~marks) Ysim <- rmh(fit) plot(Ysim, main="simulation from ppm(Y, ~marks)") fit <- ppm(Y,~polynom(x,y,2)) Ysim <- rmh(fit) plot(Ysim, main="simulation from ppm(Y, ~polynom(x,y,2))") fit <- ppm(Y,~marks+polynom(x,y,2)) Ysim <- rmh(fit) plot(Ysim, main="simulation from ppm(Y, ~marks+polynom(x,y,2))") fit <- ppm(Y,~marks*polynom(x,y,2)) Ysim <- rmh(fit) plot(Ysim, main="simulation from ppm(Y, ~marks*polynom(x,y,2))") # multitype Strauss models MS <- MultiStrauss(types = levels(Y$marks), radii=matrix(0.07, ncol=2, nrow=2)) fit <- ppm(Y, ~marks, MS) Ysim <- rmh(fit, control=list(nrep=1e3)) plot(Ysim, main="simulation from fitted Multitype Strauss") fit <- ppm(Y,~marks*polynom(x,y,2), MS) Ysim <- rmh(fit, control=list(nrep=1e3)) plot(Ysim, main="simulation from fitted inhomogeneous Multitype Strauss") ## End(Not run)