estgtp function

Bayesian MCMC estimation of parameters of generalized Thomas process

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

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

Returns

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 = "")