# rmh

##### Simulate point patterns using the Metropolis-Hastings algorithm.

Simulates point patterns from a (still small) range of point process models.

- Keywords
- spatial

##### Usage

```
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)
```

##### Arguments

- 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),`str`

- 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-rectangu - 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 shou - 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 ter
- 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 - 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 - 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`

i - 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.

##### Details

The parameter vector (or list) `par`

should be
as follows, for each of the available conditional intensity
functions:
[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object]

##### Value

- 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`

.

##### Warnings

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.

##### References

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

##### See Also

`ppp`

,
`mpl`

,
`Strauss`

,
`Softcore`

,
`StraussHard`

,
`MultiStrauss`

,
`MultiStraussHard`

##### Examples

```
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)
<testonly>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=100,nverb=50)
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=100,nverb=50)
# 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=100,nverb=50,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=100,nverb=50)
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=100,nverb=50)
# 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=100,nverb=50,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=100,
nverb=50)
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=100,nverb=50)
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=100,nverb=50)
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=100,nverb=50,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=100,nverb=50,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=100,nverb=50)
par8 <- c(1800,3,0.02,0.04)
X.dig2 <- rmh("dig2",par=par8,w=w,n.start=300,nrep=100,nverb=50)
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=100,nverb=50)
# 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=100,nverb=50)</testonly>
```

*Documentation reproduced from package spatstat, version 1.0-1, License: GPL version 2 or newer*