[R-sig-Geo] Simulating a Gibbs Marked point process

[Adrian.Baddeley at csiro.au]

As Rolf says below, the function 'rmh' in spatstat generates *simulated realisations* of a Gibbs point process.
It does not *maximise* the probability density of the Gibbs point process. 

Examples of the use of rmh are provided in the help files for rmh.default and rmh.ppm, and in the workshop notes <www.csiro.au/resources/pf16h.html>

Currently there is no command in spatstat for maximising the Gibbs point process probability density (for example by Simulated Annealing, or Steepest Ascent).

However it is easy to implement a simple form of Simulated Annealing for a specific model. If f(x) is the probability density of the model, then for each 'temperature' T > 0 you have to figure out the parameters of the modified model with probability density f_T(x) proportional to f(x)^(1/T). For example for the  stationary Strauss process with parameters beta and gamma, the modified parameters are beta_T = beta^(1/T) and gamma = gamma^(1/T). Now suppose your annealing schedule is to run at temperature T1 for n1 iterations, then temperature T2 < T1 for n2 iterations, etc. To do the first run, call 'rmh' with the appropriate parameter values for the model, with nrep=n1. Save the resulting point pattern as X1. To do the second run, call 'rmh' again with the new parameter values appropriate to temperature T2, and using the previous state as the initial condition: start=list(x.start=X1). And so on.

An important question is whether this algorithm will converge to a unique solution. There could be many point patterns which are maximal for the probability density. For example if the model is the stationary Strauss process then the limit of the probability density as T -> 0 is the density of a Hard Core process, which has many maximal configurations. If the model is non-stationary, then it may have a unique maximal configuration, but that cannot be assumed in general.

Prof Adrian Baddeley FAA
School of Earth & Environment, University of Western Australia, 35 Stirling Hwy, Crawley WA 6009, Australia
                                            - and -
CSIRO Mathematics, Informatics & Statistics, Leeuwin Centre, 65 Brockway Rd, Floreat WA 6014, Australia
skype adrian.baddeley
________________________________________
From: Rolf Turner [r.turner at auckland.ac.nz]
Sent: Wednesday, 10 April 2013 5:13 PM
To: Ferra Xu
Cc: r-sig-geo at r-project.org; Baddeley, Adrian (CMIS, Floreat)
Subject: Re: [R-sig-Geo] Simulating a Gibbs Marked point process

On 04/10/2013 02:30 PM, Ferra Xu wrote:
> Hi everybody
>
> I want to simulate a pattern of points based on maximization of Gibbs marked point process density function by MCMC methods using R. I already know that rmh function in spatstat will do that, but I couldn't find any good example code for that... I don't know how should I determine the objective function that I want to be maximized in each iteration...
I don't really understand what you want to do.  The rmh() [random
Metropolis Hastings]
function does not proceed by maximizing the density function of the
(Gibbs) process.  It
uses a Markov chain (whose "states" are point patterns), the
Metropolis-Hastings algorithm,
and the *conditional intensity* function of the process, to effect the
simulation.  The transitions
between states are governed by transition probabilities constructed so
that the *steady
state* distribution of the Markov chain is the distribution of the
desired process.  There is
no maximization of anything --- at least not in any obvious way, as far
as I can see.

What do you mean by "I couldn't find any good example code"?  There are
many examples
of code *using* the rmh() function in the online help.  If you want to
see the code of rmh()
--- well, R and all of its contributed packages are open source.
Download the spatstat source
package and have a look! (Be warned --- the heavy lifting is done in C,
and the code is pretty
subtle and obscure.  By my standards.)

     cheers,

         Rolf Turner