estgtp {binspp}

Bayesian MCMC estimation of parameters of generalized Thomas process

Description

Bayesian MCMC estimation of parameters of generalized Thomas process. The cluster size is allowed to have a variance that is greater or less than the expected value (cluster sizes are over or under dispersed).

Usage

estgtp(
  X,
  kappa0 = exp(a_kappa + ((b_kappa^2)/2)),
  omega0 = exp(a_omega + ((b_omega^2)/2)),
  lambda0 = (l_lambda + u_lambda)/2,
  theta0 = exp(a_theta + ((b_theta^2)/2)),
  skappa,
  somega,
  dlambda,
  stheta,
  smove,
  a_kappa,
  b_kappa,
  a_omega,
  b_omega,
  l_lambda,
  u_lambda,
  a_theta,
  b_theta,
  iter = 5e+05,
  plot.step = 1000,
  save.step = 1000,
  filename
)

Arguments

X

A point pattern dataset (object of class ppp) to which the model should be fitted.

kappa0

Initial value for kappa, by default it will be set as expectation of prior for kappa.

omega0

Initial value for omega, by default it will be set as expectation of prior for omega.

lambda0

Initial value for lambda, by default it will be set as expectation of prior for lambda.

theta0

Initial value for theta, by default it will be set as expectation of prior for theta.

skappa

variability of proposal for kappa: second parameter of log-normal distribution

somega

variability of proposal for omega: second parameter of log-normal distribution

dlambda

variability of proposal for lambda: half of range of uniform distribution

stheta

variability of proposal for theta: second parameter of log-normal distribution

smove

variability of proposal for moving center point: SD of normal distribution

a_kappa

First parameter of prior distribution for kappa, which is log-normal distribution.

b_kappa

Second parameter of prior distribution for kappa, which is log-normal distribution.

a_omega

First parameter of prior distribution for omega, which is log-normal distribution.

b_omega

Second parameter of prior distribution for omega, which is log-normal distribution.

l_lambda

First parameter of prior distribution for lambda, which is uniform distribution.

u_lambda

Second parameter of prior distribution for lambda, which is uniform distribution.

a_theta

First parameter of prior distribution for theta, which is log-normal distribution.

b_theta

Second parameter of prior distribution for theta, which is log-normal distribution.

iter

Number of iterations of MCMC.

plot.step

Step for the graph plotting. If the value is greater than iter parameter value, no plots will be visible.

save.step

Step for the parameters saving. The file must be specified or has to be set to larger than iter.

filename

The name of the output RDS file

Value

The output is an estimated MCMC chain of parameters, centers and connections.

Examples

library(spatstat)
kappa = 10
omega = .1
lambda= .5
theta = 10

X = rgtp(kappa, omega, lambda, theta, win = owin(c(0, 1), c(0, 1)))
plot(X$X)
plot(X$C)

a_kappa = 4
b_kappa = 1
x <- seq(0, 100, length = 100)
hx <- dlnorm(x, a_kappa, b_kappa)
plot(x, hx, type = "l", lty = 1, xlab = "x value",
     ylab = "Density", main = "Prior")

a_omega = -3
b_omega = 1
x <- seq(0, 1, length = 100)
hx <- dlnorm(x, a_omega, b_omega)
plot(x, hx, type = "l", lty = 1, xlab = "x value",
     ylab = "Density", main = "Prior")

l_lambda = -1
u_lambda = 0.99
x <- seq(-1, 1, length = 100)

hx <- dunif(x, l_lambda, u_lambda)
plot(x, hx, type = "l", lty = 1, xlab = "x value",
     ylab = "Density", main = "Prior")

a_theta = 4
b_theta = 1
x <- seq(0, 100, length = 100)
hx <- dlnorm(x, a_theta, b_theta)
plot(x, hx, type = "l", lty = 1, xlab = "x value",
     ylab = "Density", main = "Prior")

est = estgtp(X$X,
          skappa = exp(a_kappa + ((b_kappa ^ 2) / 2)) / 100,
          somega = exp(a_omega + ((b_omega ^ 2) / 2)) / 100,
          dlambda = 0.01,
          stheta = exp(a_theta + ((b_theta ^ 2) / 2)) / 100, smove = 0.1,
          a_kappa = a_kappa, b_kappa = b_kappa,
          a_omega = a_omega, b_omega = b_omega,
          l_lambda = l_lambda, u_lambda = u_lambda,
          a_theta = a_theta, b_theta = b_theta,
          iter = 50, plot.step = 50, save.step = 1e9,
          filename = "")