rmh.default(model,start,control,verbose=TRUE,...)
n.start
and x.start
are incompatible;
precisely one of them should be specifverbose
is TRUE then
warnings are printed out whenever the storage space alloted to
the underlying Fortran code, to hold the generated points, gets
increased. If it is FALSE then this process proceeds silently."ppp"
, see
ppp.object
). The returned 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, trend)
start=list(n.start,x.start,iseed)
control=list(p=p,q=q,nrep=nrep,expand,periodic,
ptypes=ptypes,fixall=fixall)
Note that only one of x.start
and x.start
appear in
in the start
list.
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:
[object Object],[object Object]
The argument start
determines the initial state of the
Metropolis-Hastings algorithm. Possible components are
[object Object],[object Object],[object Object]
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:
[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object]
Diggle, P. J. (2003) Statistical Analysis of Spatial Point Patterns (2nd ed.) Arnold, London.
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{
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. }
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
.The syntax of rmh.default
in respect of the lookup
cif
has changed from the previous release of spatstat
(versions 1.4-7 to 1.5-1).
Read the Details carefully. In particular it is now required
that the first entry of the r
component of par
be strictly positive. (This is the opposite of what was
required in the previous release, which was that this first entry
had to be 0.)
It is also now required that the entries of r
be
sorted into ascending order. (In the previous release it
was assumed that the entries of r
and h
were
in corresponding order and the two vectors were sorted
commensurately. It was decided that this is dangerous sand
unnecessary.)
Note that if you specify the lookup
pairwise interaction
function via stepfun()
the arguments x
and y
which are passed to stepfun()
are slightly
different from r
and h
: length(y)
is equal
to 1+length(x)
; the final entry of y
must be equal
to 1 --- i.e. this value is explicitly supplied by the user rather
than getting tacked on internally.
The step function returned by stepfun()
must be right
continuous (this is the default behaviour of stepfun()
)
otherwise an error is given.
}
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 x.start$window
. Any specified
value for expand
is simply ignored.
The presence of both a component w
of model
and a
non-null value for x.start$window
makes sense ONLY if w
is contained in x.start$window
. However no checking
is done for this.
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.
No checking is done to verify non-negativity of any specified trend or trends. }
rmh
,
rmh.ppm
,
ppp
,
ppm
,
Strauss
,
Softcore
,
StraussHard
,
MultiStrauss
,
MultiStraussHard
,
DiggleGratton
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
.)
Note that the lookup
conditional intensity function
permits the simulation (in theory, to any desired degree
of approximation) of any pairwise interaction process for
which the interaction depends only on the distance between
the pair of points.
}
# Lookup (interaction function h_2 from page 76, Diggle (2003)):
r <- seq(from=0,to=0.2,length=101)[-1] # Drop 0.
h <- 20*(r-0.05)
h[r<0.05] 0="" <-="" h[r="">0.10] <- 1
mod17 <- list(cif="lookup",par=list(beta=4000,h=h,r=r),w=c(0,1,0,1))
X.lookup <- rmh(model=mod17,start=list(n.start=100),
control=list(nrep=nr,nverb=nv))
# Strauss with trend
tr <- function(x,y){x <- x/250; y <- y/250;
exp((6*x + 5*y - 18*x^2 + 12*x*y - 9*y^2)/6)
}
beta <- 0.3
gmma <- 0.5
r <- 45
mod17 <- list(cif="strauss",par=c(beta=beta,gamma=gmma,r=r),w=c(0,250,0,250),
trend=tr3,tmax=1.5)
X1.strauss.trend <- rmh(model=mod17,start=list(n.start=90),
control=list(nrep=nr,nverb=nv))
[object Object],[object Object]