rmh                 package:spatstat                 R Documentation

_S_i_m_u_l_a_t_e _p_o_i_n_t _p_a_t_t_e_r_n_s _u_s_i_n_g _t_h_e _M_e_t_r_o_p_o_l_i_s-_H_a_s_t_i_n_g_s _a_l_g_o_r_i_t_h_m.

_D_e_s_c_r_i_p_t_i_o_n:

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

_U_s_a_g_e:

     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)

_A_r_g_u_m_e_n_t_s:

     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), `straush' (Strauss process
          with hardcore), `sftcr' (Softcore process), `straussm'
          (multitype Strauss process), `straushm' (multitype Strauss
          process with hardcore), `dig1', `dig2' (see the 2nd of the
          References), and `geyer' (see the 3rd of the References).

          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 or ratfor code for the new
          conditional intensity function, the user would have to modify
          the code in the files `cif.r' and `rmh.R' appropriately. 
          (And recreate/recompile the dynamically loadable shared
          object `spatstat.so'.) 

     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-rectangular windows we cannot (as yet) generated a
          pattern from a process which exists only within that window. 
          The domain of the theoretical process must (at present) be at
          least as large as the ``enclosing box'' of the window. 

  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 should be improved if
          `ptypes' is close to the relative frequencies of the types
          which will result from the simulation. 

    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 terms x, y, x^2, xy, y^2, x^3, x^2y, xy^2, y^3, ..., y^n
          where n is the degree of the trend.  (Thus the possible
          lengths of a sequence of coefficients are 0, 2, 5, 9, 14,
          ...).

          If the process is unmarked `tpar' must be a vector consisting
          solely of the coefficients of the trend polynomial in the
          correct order.  The degree of the polynomial will be
          calculated from the length of `tpar'.

          If the process is multitype and if `tpar' is given as a
          vector, then `ntypes' and and the degrees of each of the
          polynomial should comprise the first 1+`ntypes' entries of
          `tpar'.  After that, the coefficients corresponding to each
          mark should be catenated successively.  Note that some of the
          degrees are allowed to be 0.

          If `tpar' is given as a list, the i-th component of this list
          should consist of the coefficients for the trend
          corresponding to the i-th type.  If there is to be no trend
          for the i-th type, then the i-th component of the list should
          be NULL (but must be present in the list as the NULL object.) 

 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 `n.start'
          gets multiplied by the ratio of the area of the enclosing box
          for the window to the area of the window, and then by the
          factor `expand'.  Then that many points are uniformly
          generated in the expanded window; see below.)  The value of
          `n.start' should be roughly equal to (an educated guess at)
          the expected number of points which will be generated inside
          the window. 

  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 1, then we are
          simulating the latter (unless `periodic' in `TRUE'; see
          below).  The larger `expand' is, the better we approximate
          the former.  Note that any value of `expand' smaller than 1
          is treated as if it were 1. The area of the expanded window
          is equal to `expand' times the area of the enclosing box;
          width and height are stretched proportionately.  Points are
          generated by the Metropolis-Hastings algorithm in the
          expanded window, and then ``clipped'' down to the original
          window `w' when the algorithm has finished. The argument
          `expand' defaults to 2 when `periodic' is `FALSE' and to 1
          when periodic is `TRUE'.  (Trying to set `expand' greater
          than 1 when periodic is `TRUE' generates an error.) 

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' is not rectangular, an error is given. 

    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. 

_D_e_t_a_i_l_s:

     The parameter vector (or list) `par' should be as follows, for
     each of the available conditional intensity functions:

     _s_t_r_a_u_s_s: (Strauss process.) A vector with components beta,gamma,r
          which are respectively the ``base'' intensity, the pair-wise
          interaction parameter and the interaction radius.  Note that
          gamma must be less than or equal to 1.

     _s_t_r_a_u_s_h: (Strauss process with hardcore.) A vector with entries
          beta,gamma,r,r_hc where beta, gamma, and r are as for the
          Strauss process, and r_hc is the hardcore radius.  Of course
          r_hc must be less than r.

     _s_f_t_c_r: (Softcore process.) A vector with components
          beta,sigma,kappa.  Again beta is a ``base'' intensity.  The
          pairwise interaction between two points u != v is

                      -(sigma/||u-v||)^(2/kappa)

          Note that it is necessary that 0 < kappa <1.

     _s_t_r_a_u_s_s_m: (Multitype Strauss process.) Here `par' is best given as
          a list with components

          _b_e_t_a: A vector of ``base'' intensities, one for each possible
               type.

          _g_a_m_m_a: A symmetric matrix of interaction parameters, with
               gamma_ij pertaining to the interaction between type i
               and type j.

          _r: A symmetric matrix of interaction radii, with r_ij
               pertaining to the interaction between type i and type j.

          If `par' is to be given as a vector, then this vector should
          consist of the vector `beta', catenated with the upper
          triangle of the matrix `gamma' strung out in row order, and
          finally catenated with the upper triangle of the matrix `r'
          likewise strung out in row order.

     _s_t_r_a_u_s_h_m: (Multitype Strauss process with hardcore.) The
          parameters `par' are much as for `straussm' except that there
          is an extra component `rhc' which is the matrix of hardcore
          radii.  Note that in the vector form, the hardcore radii must
          be catenated after the interaction radii.

     _d_i_g_1: (See the 2nd of the References.) Process with pairwise
          interaction function

                     e(t) = sin^2((pi t)/(2 rho))

          for t < rho, and equal to 1 for t >= rho.  The parameters are
          beta and rho.

     _d_i_g_2: (See the 2nd of the References.) Process with pairwise
          interaction function e(t) equal to 0 for t < delta, equal to

                    ((t-delta)/(rho-delta))^kappa

          for delta <= t < rho, and equal to 1 for delta >= rho.  Note
          that here we use the symbol kappa where Diggle, Gates, and
          Stibbard use beta since we reserve the symbol beta for an
          intensity parameter.

          The complete parameter set is beta, kappa, delta and rho.

     _g_e_y_e_r (See the 3rd of the References.) Geyer's ``saturation''
          point process model.  This model is ``like a Strauss model,
          but with an upper bound to the number of r-close neighbors of
          any point.''

          More explicitly, a saturation point process with interaction
          radius r, saturation threshold s, and  parameters beta and
          gamma, is the point process in which each point x_i in the
          pattern X contributes a factor

                      beta gamma^min(s,t(x_i,X))

          to the probability density of the point pattern, where
          t(x_i,X) denotes the number of ``r-close neighbours'' of x_i
          in the pattern X.

_V_a_l_u_e:

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

_W_a_r_n_i_n_g_s:

     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.

_A_u_t_h_o_r(_s):

     Adrian Baddeley adrian@maths.uwa.edu.au <URL:
     http://www.maths.uwa.edu.au/~adrian/> and Rolf Turner
     rolf@math.unb.ca <URL: http://www.math.unb.ca/~rolf>

_R_e_f_e_r_e_n_c_e_s:

     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

_S_e_e _A_l_s_o:

     `ppp', `mpl', `Strauss', `Softcore', `StraussHard',
     `MultiStrauss', `MultiStraussHard'

_E_x_a_m_p_l_e_s:

        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)



