Nothing
#'
#' clustersiminfo.R
#'
#' Copyright (C) Adrian Baddeley and Ya-Mei Chang 2022-2023
#' GNU Public Licence >= 2
#'
#' -----------------------------------------------------------------
#'
#' Table of additional information for cluster models
#' for use in Brix-Kendall and Baddeley-Chang simulation algorithms
#'
#' D = bounding disc
#' rD = radius of D (numeric > 0)
#' r = distance from parent point to centre of disc (numeric vector)
#' mod = list containing all model parameters:
#' $par = original (native) parameters of cluster model e.g. (kappa, sigma2)
#' $mu = mean number of offspring per parent in stationary process
#' $margs = list of shape parameters
#'
#' Entries:
#'
#' (A) cluster model
#'
#' roffspring = function(n, mod)
#' generate x, y coordinates of n offspring of parent at origin
#'
#' (A) Brix-Kendall dominating process
#'
#' hdom = function(r, mod, rD)
#' Value inside D of the dominating offspring density, for a parent at distance r from centre of D
#'
#' Eplus = function(r, mod, rD)
#' Integral over D of dominating kernel, for a parent at distance r from centre of D
#' Eplus(r) = pi * rD^2 * mu * hdom(r)
#'
#' rhoplus = function(r, mod, rD)
#' Intensity of Brix-Kendall dominating parent process, given distance to origin
#' rhoplus(r) = kappa * (1 - exp(-Eplus(r)))
#'
#' Mplus = function(r, mod, rD, ...)
#' Radial cumulative integral of intensity of Brix-Kendall dominating intensity
#' Mplus(r) = \int_0^r 2 pi t rhoplus(t) dt
#'
#' MplusInf = function(mod, rD)
#' Total integral Mplus(infty)
#' Expected total number of parents in the dominating process
#'
#' invMplus = function(v, mod, rD, Minfty)
#' Inverse function of Mplus (IF KNOWN)
#'
#' (B) Baddeley-Chang super-dominating process
#'
#' rhoplusplus = function(r, mod, rD)
#' Intensity of superdominating process
#' rhoplusplus(r) = kappa * Eplus(r)
#'
#' Mplusplus = function(r, mod, rD, ...)
#' Radial cumulative integral of superdominating intensity
#' Mplusplus(r) = \int_0^r 2 pi t rhoplusplus(t) dt
#' = kappa \int_0^r 2 pi t Eplus(t) dt
#'
#' MplusplusInf = function(mod, rD)
#' Total integral Mplusplus(infty)
#' Expected total number of parents in the superdominating process
#'
#' invMplusplus = function(v, mod, rD, Minfty)
#' Inverse function of Mplusplus (IF KNOWN)
#'
#' inflate = function(mod, rD)
#' Rule for determining optimal inflated radius rE
#' according to Baddeley and Chang (2023) section 6.6
#' >>>>>>>>>>>>>>>> THOMAS PROCESS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
.Thomas.roffspring <- function(n, mod) {
sigma <- sqrt(mod$par[["sigma2"]])
list(x=rnorm(n, sd=sigma), y=rnorm(n, sd=sigma))
}
.Thomas.hdom <- function(r, mod, rD) {
inv2sig2 <- 1/(2 * mod$par[["sigma2"]])
B <- inv2sig2/pi
z <- numeric(length(r))
above <- (r > rD)
z[!above] <- B
z[above] <- B * exp(- inv2sig2 * (r[above] - rD)^2)
return(z)
}
.Thomas.Eplus <- function(r, mod, rD) {
mu <- mod$mu
inv2sig2 <- 1/(2 * mod$par[["sigma2"]])
A <- mu * rD^2 * inv2sig2
z <- numeric(length(r))
above <- (r > rD)
z[!above] <- A
z[above] <- A * exp(- inv2sig2 * (r[above]-rD)^2)
return(z)
}
.Thomas.rhoplus <- function(r, mod, rD) {
kappa <- mod$par[["kappa"]]
z <- numeric(length(r))
above <- (r > rD)
z[!above] <- kappa * (1 - exp(-.Thomas.Eplus(0, mod, rD)))
z[above] <- kappa * (1 - exp(-.Thomas.Eplus(r[above], mod, rD)))
return(z)
}
.Thomas.MplusIntegrand <- function(r, mod, rD) {
2 * pi * r * .Thomas.rhoplus(r, mod, rD)
}
.Thomas.Mplus <- function(r, mod, rD, ..., method="q") {
## integrand is linear on [0, rD]
z <- pi * pmin(r, rD)^2 * .Thomas.rhoplus(0, mod, rD)
high <- (r > rD)
if(any(high))
z[high] <- z[high] + indefinteg(.Thomas.MplusIntegrand,
r[high],
lower=rD,
mod=mod, rD=rD, method=method)
return(z)
}
.Thomas.MplusInf <- function(mod, rD) { .Thomas.Mplus(Inf, mod, rD) }
## inverse function requires numerical root-finding
.Thomas.rhoplusplus <- function(r, mod, rD) {
mod$par[["kappa"]] * .Thomas.Eplus(r, mod, rD)
}
.Thomas.Mplusplus <- function(r, mod, rD, ...) {
mu <- mod$mu
kappa <- mod$par[["kappa"]]
sigma2 <- mod$par[["sigma2"]]
inv2sig2 <- 1/(2 * sigma2)
z <- pi * pmin(r, rD)^2
high <- (r > rD)
z[high] <- z[high] +
2 * pi * sigma2 * (1 - exp(-inv2sig2 * (r[high]-rD)^2)) +
2 * pi * rD * sqrt(2 * pi * sigma2) * (
pnorm(r[high]-rD, sd=sqrt(sigma2)) - 1/2
)
z <- kappa * mu * rD^2 * inv2sig2 * z
return(z)
}
.Thomas.MplusplusInf <- function(mod, rD) {
mu <- mod$mu
kappa <- mod$par[["kappa"]]
sigma2 <- mod$par[["sigma2"]]
inv2sig2 <- 1/(2 * sigma2)
kappa * mu * pi * rD^2 * (
rD^2 * inv2sig2 +
1 +
rD * sqrt(2 * pi * sigma2) * inv2sig2
)
}
## inverse function requires numerical root-finding
.Thomas.inflate <- function(mod, rD) {
mu <- mod$mu
sigma2 <- mod$par[["sigma2"]]
a <- if(mu == 0) 1 else (1 + (1-exp(-mu))/mu)
b <- (rD^2)/(2*a*sigma2)
if(b <= 1) return(rD)
delta <- 2 * sqrt(sigma2) * sqrt(log(b)/2)
return(rD + delta)
}
#' >>>>>>>>>>>>>>>> MATERN CLUSTER PROCESS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
.MatClust.roffspring <- function(n, mod) {
R <- mod$par[["R"]]
rad <- R * sqrt(runif(n))
theta <- runif(n, max=2*pi)
list(x = rad * cos(theta), y = rad * sin(theta))
}
.MatClust.hdom <- function(r, mod, rD) {
R <- mod$par[["R"]]
(r - rD < R)/(pi * R^2)
}
.MatClust.Eplus <- function(r, mod, rD) {
R <- mod$par[["R"]]
mu <- mod$mu
mu * (r - rD < R) * (rD/R)^2
}
.MatClust.rhoplus <- function(r, mod, rD) {
mod$par[["kappa"]] * (1 - exp(-.MatClust.Eplus(r, mod, rD)))
}
.MatClust.Mplus <- function(r, mod, rD, ...) {
R <- mod$par[["R"]]
kappa <- mod$par[["kappa"]]
mu <- mod$mu
pi * pmin(r, R+rD)^2 * kappa * ( 1 - exp(-mu * (rD/R)^2))
}
.MatClust.MplusInf <- function(mod, rD) {
R <- mod$par[["R"]]
kappa <- mod$par[["kappa"]]
mu <- mod$mu
pi * (R+rD)^2 * kappa * ( 1 - exp(-mu * (rD/R)^2))
}
.MatClust.inverseMplus <- function(v, mod, rD, Minfty) {
R <- mod$par[["R"]]
kappa <- mod$par[["kappa"]]
mu <- mod$mu
B <- pi * kappa * (1 - exp(- mu * rD^2/R^2))
sqrt(v/B)
}
.MatClust.rhoplusplus <- function(r, mod, rD) {
R <- mod$par[["R"]]
kappa <- mod$par[["kappa"]]
mu <- mod$mu
kappa * mu * (r < (R+rD)) * (rD/R)^2
}
.MatClust.Mplusplus <- function(r, mod, rD, ...) {
R <- mod$par[["R"]]
kappa <- mod$par[["kappa"]]
mu <- mod$mu
kappa * mu * pi * (pmin(r, R+rD) * rD/R)^2
}
.MatClust.MplusplusInf <- function(mod, rD) {
R <- mod$par[["R"]]
kappa <- mod$par[["kappa"]]
mu <- mod$mu
kappa * mu * pi * ((R+rD) * rD/R)^2
}
.MatClust.inverseMplusplus <- function(v, mod, rD, Minfty) {
R <- mod$par[["R"]]
kappa <- mod$par[["kappa"]]
mu <- mod$mu
(R/rD) * sqrt(v/(pi * kappa * mu))
}
.MatClust.inflate <- function(mod, rD) {
R <- mod$par[["R"]]
rE <- if(R < rD) (rD + R) else rD
return(rE)
}
#' >>>>>>>>>>>>>>>> CAUCHY CLUSTER PROCESS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
#' idiosyncratic parameter eta2 = 4 * scale^2
#' scale <- sqrt(mod$par[["eta2"]])/2
.Cauchy.roffspring <- function(n, mod) {
rate <- mod$par[["eta2"]]/8
b <- 1/sqrt(rgamma(n, shape=1/2, rate=rate))
list(x = b * rnorm(n), y = b * rnorm(n))
}
.Cauchy.hdom <- function(r, mod, rD) {
scale2 <- mod$par[["eta2"]]/4
z <- numeric(length(r))
above <- (r > rD)
z[!above] <- 1/(2*pi*scale2)
z[above] <- (1/(2*pi*scale2)) * (1 + ((r[above]-rD)^2)/scale2)^(-3/2)
return(z)
}
.Cauchy.Eplus <- function(r, mod, rD) {
scale2 <- mod$par[["eta2"]]/4
mu <- mod$mu
z <- numeric(length(r))
above <- (r > rD)
A <- mu * (rD^2/(2*scale2))
z[!above] <- A
z[above] <- A * (1 + ((r[above]-rD)^2)/scale2)^(-3/2)
return(z)
}
.Cauchy.rhoplus <- function(r, mod, rD) {
kappa <- mod$par[["kappa"]]
z <- numeric(length(r))
above <- (r > rD)
z[!above] <- kappa * (1 - exp(-.Cauchy.Eplus(0, mod, rD)))
z[above] <- kappa * (1 - exp(-.Cauchy.Eplus(r[above], mod, rD)))
return(z)
}
.Cauchy.MplusIntegrand <- function(r, mod, rD) {
2 * pi * r * .Cauchy.rhoplus(r, mod, rD)
}
.Cauchy.Mplus <- function(r, mod, rD, ..., method="q") {
## integrand is linear on [0, rD]
z <- pi * pmin(r, rD)^2 * .Cauchy.rhoplus(0, mod, rD)
high <- (r > rD)
if(any(high))
z[high] <- z[high] + indefinteg(.Cauchy.MplusIntegrand,
r[high],
lower=rD,
mod=mod, rD=rD, method=method)
return(z)
}
.Cauchy.MplusInf <- function(mod, rD) { .Cauchy.Mplus(Inf, mod, rD) }
.Cauchy.rhoplusplus <- function(r, mod, rD) {
mod$par[["kappa"]] * .Cauchy.Eplus(r, mod, rD)
}
.Cauchy.Mplusplus <- function(r, mod, rD, ...) {
mu <- mod$mu
kappa <- mod$par[["kappa"]]
lambda <- kappa * mu
scale2 <- mod$par[["eta2"]]/4
distD <- pmax(0, r - rD)
lambda * pi * rD^2 * ifelse(r <= rD,
r^2/(2*scale2),
rD^2/(2*scale2) + 1 + (rD * distD/scale2 - 1)/sqrt(1+distD^2/scale2))
}
.Cauchy.MplusplusInf <- function(mod, rD) {
mu <- mod$mu
kappa <- mod$par[["kappa"]]
lambda <- kappa * mu
scale2 <- mod$par[["eta2"]]/4
rdons2 <- rD^2/scale2
lambda * pi * rD^2 * (1 + sqrt(rdons2) + rdons2/2)
}
.Cauchy.inflate <- function(mod, rD) {
mu <- mod$mu
scale2 <- mod$par[["eta2"]]/4
a <- if(mu == 0) 1 else (1 + (1-exp(-mu))/mu)
b <- (rD^2)/(2*a*scale2)
if(b <= 1) return(rD)
delta <- sqrt(scale2) * (b^(2/3) - 1)
return(rD + delta)
}
#' ........................................................................
.spatstat.clustersim.InfoTable <- list(
Thomas = list(
iscompact = FALSE,
roffspring = .Thomas.roffspring,
hdom = .Thomas.hdom,
Eplus = .Thomas.Eplus,
rhoplus = .Thomas.rhoplus,
Mplus = .Thomas.Mplus,
MplusInf = .Thomas.MplusInf,
invMplus = NULL,
rhoplusplus = .Thomas.rhoplusplus,
Mplusplus = .Thomas.Mplusplus,
MplusplusInf = .Thomas.MplusplusInf,
invMplusplus = NULL,
inflate = .Thomas.inflate
),
MatClust = list(
iscompact = TRUE,
roffspring = .MatClust.roffspring,
hdom = .MatClust.hdom,
Eplus = .MatClust.Eplus,
rhoplus = .MatClust.rhoplus,
Mplus = .MatClust.Mplus,
MplusInf = .MatClust.MplusInf,
invMplus = .MatClust.inverseMplus,
rhoplusplus = .MatClust.rhoplusplus,
Mplusplus = .MatClust.Mplusplus,
MplusplusInf = .MatClust.MplusplusInf,
invMplusplus = .MatClust.inverseMplusplus,
inflate = .MatClust.inflate
),
Cauchy = list(
iscompact = FALSE,
roffspring = .Cauchy.roffspring,
hdom = .Cauchy.hdom,
Eplus = .Cauchy.Eplus,
rhoplus = .Cauchy.rhoplus,
Mplus = .Cauchy.Mplus,
MplusInf = .Cauchy.MplusInf,
invMplus = NULL,
rhoplusplus = .Cauchy.rhoplusplus,
Mplusplus = .Cauchy.Mplusplus,
MplusplusInf = .Cauchy.MplusplusInf,
invMplusplus = NULL,
inflate = .Cauchy.inflate
)
)
spatstatClusterSimModelMatch <- function(name, verbose=TRUE) {
if(!is.character(name) || length(name) != 1)
stop("Argument must be a single character string", call.=FALSE)
TheTable <- .spatstat.clustersim.InfoTable
nama2 <- names(TheTable)
mat <- pmatch(name, nama2)
if(!is.null(mat))
return(nama2[mat])
if(verbose)
warning(sQuote(name), "is not supported;",
"available options are", commasep(sQuote(nama2)))
return(NULL)
}
spatstatClusterSimInfo <- function(name) {
known <- spatstatClusterSimModelMatch(name, FALSE)
if(is.null(known)) return(NULL)
return(.spatstat.clustersim.InfoTable[[known]])
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.