mple.cppm {NScluster}

MPLE of Neyman-Scott Cluster Point Process Models and Their Extensions

Description

MPLE of the five cluster point process models.

Usage

mple.cppm(model = "Thomas", xy.points, pars = NULL, eps = 0.001, uplimit = 0.3,
          skip = 1)

## S3 method for class 'mple'
coef(object, ...)
## S3 method for class 'mple'
summary(object, ...)

Arguments

model

a character string indicating each cluster point process model: "Thomas", "IP", "TypeA", "TypeB", and "TypeC".

xy.points

a matrix containing the coordinates (x,y) of points in W=[0,1]\times[0,1].

pars

a named vector containing a given initial guess of each parameter. If NULL, any suitable parameters are used. See Details' in sim.cppm for the parameters of each model.

eps

the sufficiently small number to implement the optimization procedure for the log-Palm likelihood function. The procedure is iterated at most 1000 times until the process2$stderr becomes smaller than the eps.

uplimit

upper limit in place of \infty of the integral in the probability distribution function relative to the random distance between two descendant points within the same cluster. The uplimit is valid for "IP" and "TypeA".

skip

the variable enables one to obtain speedily the initial MPLEs, but rough approximation. The skip calculates the Palm intensity function of the log-Palm likelihood function for every skip-th r_{ij} in the ordered distances of the pairs i and j. The skip is valid for "IP" and "TypeA".

object

an object of class "mple".

...

ignored.

Details

• "Thomas" (Thomas model)

  The Palm intensity function is given as follows:

  For all r \ge 0,

  \lambda_{\bm{o}}(r) = \mu\nu + \frac{\nu}{4\pi \sigma^2} \exp \left( -\frac{r^2}{4 \sigma^2} \right).

  The log-Palm likelihood function is given by

  \log L(\mu,\nu,\sigma) = \sum_{\{i,j; i<j, r_{ij} \le 1/2\}} \log \nu \left\{ \mu + \frac{1}{4 \pi \sigma^2} \exp \left( -\frac{{r_{ij}}^2}{4 \sigma^2} \right) \right\}

  - N(W)\nu \left\{ \frac{\pi \mu}{4} + 1 - \exp \left( -\frac{1}{16 \sigma^2} \right) \right\}.

• "TypeB" (Type B model)

  The Palm intensity function is given as follows:

  For all r \ge 0,

  \lambda_{\bm{o}}(r) = \lambda + \frac{\nu}{4 \pi} \left\{ \frac{a}{{\sigma_1}^2} \exp \left( -\frac{r^2}{4{\sigma_1}^2} \right)+ \frac{(1-a)}{{\sigma_2}^2} \exp \left( -\frac{r^2}{4{\sigma_2}^2} \right) \right\},

  where \lambda = \nu(\mu_1+\mu_2) and a = \mu_1/(\mu_1+\mu_2) are the total intensity and the ratio of the intensity of the parent points of the smaller cluster to the total one, respectively.

  The log-Palm likelihood function is given by

  \log L(\lambda, \alpha, \beta, \sigma_1, \sigma_2)

  =\sum_{\{i,j; i<j, r_{ij} \le 1/2\}} \log \left[ \lambda + \frac{1}{4 \pi} \left\{ \frac{\alpha}{{\sigma_1}^2} \exp \left( -\frac{{r_{ij}}^2}{4{\sigma_1}^2} \right) + \frac{\beta}{{\sigma_2}^2} \exp \left( -\frac{{r_{ij}}^2}{4{\sigma_2}^2} \right) \right\} \right]

  - N(W) \left[ \frac{\pi \lambda}{4} + \alpha \left\{ 1 - \exp \left( -\frac{1}{16{\sigma_1}^2} \right) \right\} + \beta \left\{ 1- \exp \left( -\frac{1}{16{\sigma_2}^2} \right) \right\} \right],

  where \alpha = a\nu and \beta = (1-a)\nu.

• "TypeC" (Type C model)

  The Palm intensity function is given as follows:

  For all r \ge 0,

  \lambda_{\bm{o}}(r) = \lambda + \frac{1}{4 \pi} \left\{ \frac{a\nu_1}{{\sigma_1}^2} \exp \left( -\frac{r^2}{4{\sigma_1}^2} \right) + \frac{(1-a)\nu_2}{{\sigma_2}^2} \exp \left( -\frac{r^2}{4{\sigma_2}^2} \right) \right\},

  where \lambda = \mu_1\nu_1 + \mu_2\nu_2 and a = \mu_1\nu_1/\lambda are the total intensity and the ratio of the intensity of the smaller cluster to the total one, respectively.

  The log-Palm likelihood function is given by

  \log L(\lambda, \alpha, \beta, \sigma_1, \sigma_2)

  = \sum_{\{i,j; i<j, r_{ij} \le 1/2\}} \log \left[ \lambda + \frac{1} {4 \pi} \left\{ \frac{\alpha}{{\sigma_1}^2} \exp \left( -\frac{{r_{ij}}^2}{4{\sigma_1}^2} \right) + \frac{\beta}{{\sigma_2}^2} \exp \left( -\frac{{r_{ij}}^2}{4{\sigma_2}^2} \right) \right\} \right]

  -N(W) \left[ \frac{\pi\lambda}{4} + \alpha \left\{ 1 - \exp \left( -\frac{1}{16{\sigma_1}^2} \right) \right\} + \beta \left\{ 1- \exp \left( -\frac{1}{16{\sigma_2}^2} \right) \right\} \right],

  where \alpha = a\nu_1 and \beta = (1-a)\nu_2.

For the inverse-power model and the Type A models, we need to take the alternative form without explicit representation of the Palm intensity function. See the second reference below for details.

Value

mple.cppm returns an object of class "mple" containing the following main components:

mple

MPLE (maximum Palm likelihood estimate).

log.mpl

the log maximum Palm likelihood.

aic

AIC.

process1

a list with following components.

cflg

1 (="update") or -1 (="testfn"), where "update" indicates that -log L value has attained the minimum so far, otherwise not.

logl

the minimized -log L in the process to minimize the negative log-Palm likelihood function.

mples

corresponding MPLEs.

process2

a list with following components.

logl.simplex

the minimized -log L by the simplex method.

stderr

standard error.

pa.normal

the normalized variables corresponding to the MPLEs relative to the initial estimates.

There are other methods plot.mple and print.mple for this class.

References

Tanaka, U., Ogata, Y. and Katsura, K. (2008) Simulation and estimation of the Neyman-Scott type spatial cluster models. Computer Science Monographs 34, 1-44. The Institute of Statistical Mathematics, Tokyo. https://www.ism.ac.jp/editsec/csm/.

Tanaka, U., Ogata, Y. and Stoyan, D. (2008) Parameter estimation and model selection for Neyman-Scott point processes. Biometrical Journal 50, 43-57.

Examples

## Not run: 
# The computation of MPLEs takes a long CPU time in the minimization procedure,
# especially for the Inverse-power type and the Type A models.

### Thomas Model
 # simulation
 pars <- c(mu = 50.0, nu = 30.0, sigma = 0.03)
 t.sim <- sim.cppm("Thomas", pars, seed = 117)
 ## estimation
 init.pars <- c(mu = 40.0, nu = 40.0, sigma = 0.05)
 t.mple <- mple.cppm("Thomas", t.sim$offspring$xy, init.pars)
 coef(t.mple)

### Inverse-Power Type Model
 # simulation
 pars <- c(mu = 50.0, nu = 30.0, p = 1.5, c = 0.005)
 ip.sim <- sim.cppm("IP", pars, seed = 353)
 ## estimation
 init.pars <- c(mu = 55.0, nu = 35.0, p = 1.0, c = 0.01)
 ip.mple <- mple.cppm("IP", ip.sim$offspring$xy, init.pars, skip = 100)
 coef(ip.mple)

### Type A Model
 # simulation
 pars <- c(mu = 50.0, nu = 30.0, a = 0.3, sigma1 = 0.005, sigma2 = 0.1)
 a.sim <- sim.cppm("TypeA", pars, seed = 575)
 ## estimation
 init.pars <- c(mu = 60.0, nu = 40.0, a = 0.5, sigma1 = 0.01, sigma2 = 0.1)
 a.mple <- mple.cppm("TypeA", a.sim$offspring$xy, init.pars, skip = 100)
 coef(a.mple)

### Type B Model
 # simulation
 pars <- c(mu1 = 10.0, mu2 = 40.0, nu = 30.0, sigma1 = 0.01, sigma2 = 0.03)
 b.sim <- sim.cppm("TypeB", pars, seed = 257)
 ## estimation
 init.pars <- c(mu1 = 20.0, mu2 = 30.0, nu = 30.0, sigma1 = 0.02, sigma2 = 0.02)
 b.mple <- mple.cppm("TypeB", b.sim$offspring$xy, init.pars)
 coef(b.mple)

### Type C Model
 # simulation
 pars <- c(mu1 = 5.0, mu2 = 9.0, nu1 = 30.0, nu2 = 150.0,
           sigma1 = 0.01, sigma2 = 0.05)
 c.sim <- sim.cppm("TypeC", pars, seed = 555)
 ## estimation
 init.pars <- c(mu1 = 10.0, mu2 = 10.0, nu1 = 30.0, nu2 = 120.0,
                sigma1 = 0.03, sigma2 = 0.03)
 c.mple <- mple.cppm("TypeC", c.sim$offspring$xy, init.pars)
 coef(c.mple)

## End(Not run)