Simulate point patterns using the Metropolis-Hastings algorithm.


Generic function for running the Metropolis-Hastings algorithm to produce simulated realisations of a point process model.


rmh(model, ...)



The point process model to be simulated.


Further arguments controlling the simulation.


The Metropolis-Hastings algorithm can be used to generate simulated realisations from a wide range of spatial point processes. For caveats, see below.

The function rmh is generic; it has methods rmh.ppm (for objects of class "ppm") and rmh.default (the default). The actual implementation of the Metropolis-Hastings algorithm is contained in rmh.default. For details of its use, see rmh.ppm or rmh.default.

[If the model is a Poisson process, then Metropolis-Hastings is not used; the Poisson model is generated directly using rpoispp or rmpoispp.]

In brief, the Metropolis-Hastings algorithm is a Markov Chain, whose states are spatial point patterns, and whose limiting distribution is the desired point process. After running the algorithm for a very large number of iterations, we may regard the state of the algorithm as a realisation from the desired point process.

However, there are difficulties in deciding whether the algorithm has run for “long enough”. The convergence of the algorithm may indeed be extremely slow. No guarantees of convergence are given!

While it is fashionable to decry the Metropolis-Hastings algorithm for its poor convergence and other properties, it has the advantage of being easy to implement for a wide range of models.


A point pattern, in the form of an object of class "ppp". See rmh.default for details.


As of version 1.22-1 of spatstat a subtle change was made to rmh.default(). We had noticed that the results produced were sometimes not “scalable” in that two models, differing in effect only by the units in which distances are measured and starting from the same seed, gave different results. This was traced to an idiosyncracy of floating point arithmetic. The code of rmh.default() has been changed so that the results produced by rmh are now scalable. The downside of this is that code which users previously ran may now give results which are different from what they formerly were.

In order to recover former behaviour (so that previous results can be reproduced) set spatstat.options(scalable=FALSE). See the last example in the help for rmh.default.



and \rolf

See Also



    # See examples in rmh.default and rmh.ppm

Want to suggest features or report bugs for Use the GitHub issue tracker. Vote for new features on Trello.

comments powered by Disqus