linearK.R
#
# linearK
#
# $Revision: 1.40 $ $Date: 2016/07/18 03:48:10 $
#
# K function for point pattern on linear network
#
#
linearK <- function(X, r=NULL, ..., correction="Ang") {
stopifnot(inherits(X, "lpp"))
correction <- pickoption("correction", correction,
c(none="none",
Ang="Ang",
best="Ang"),
multi=FALSE)
# extract info about pattern
np <- npoints(X)
lengthL <- volume(domain(X))
# compute K
denom <- np * (np - 1)/lengthL
K <- linearKengine(X, r=r, denom=denom, correction=correction, ...)
# set appropriate y axis label
switch(correction,
Ang = {
ylab <- quote(K[L](r))
fname <- c("K", "L")
},
none = {
ylab <- quote(K[net](r))
fname <- c("K", "net")
})
K <- rebadge.fv(K, new.ylab=ylab, new.fname=fname)
return(K)
}
linearKinhom <- function(X, lambda=NULL, r=NULL, ...,
correction="Ang", normalise=TRUE, normpower=1) {
stopifnot(inherits(X, "lpp"))
correction <- pickoption("correction", correction,
c(none="none",
Ang="Ang",
best="Ang"),
multi=FALSE)
if(is.null(lambda))
linearK(X, r=r, ..., correction=correction)
if(normalise) {
check.1.real(normpower)
stopifnot(normpower >= 1)
}
# extract info about pattern
lengthL <- volume(domain(X))
#
lambdaX <- getlambda.lpp(lambda, X, ...)
#
invlam <- 1/lambdaX
invlam2 <- outer(invlam, invlam, "*")
denom <- if(!normalise) lengthL else
if(normpower == 1) sum(invlam) else
lengthL * (sum(invlam)/lengthL)^normpower
K <- linearKengine(X, ...,
r=r, reweight=invlam2, denom=denom, correction=correction)
# set appropriate y axis label
switch(correction,
Ang = {
ylab <- quote(K[L, inhom](r))
fname <- c("K", "list(L, inhom)")
},
none = {
ylab <- quote(K[net, inhom](r))
fname <- c("K", "list(net, inhom)")
})
K <- rebadge.fv(K, new.fname=fname, new.ylab=ylab)
attr(K, "dangerous") <- attr(lambdaX, "dangerous")
return(K)
}
getlambda.lpp <- function(lambda, X, ..., update=TRUE) {
lambdaname <- deparse(substitute(lambda))
XX <- as.ppp(X)
missup <- missing(update)
if(update && (is.lppm(lambda) || is.ppm(lambda))) {
danger <- FALSE
lambda <- if(is.lppm(lambda)) update(lambda, X) else update(lambda, XX)
if(missup)
warn.once("lin.inhom.update",
"The behaviour of linearKinhom and similar functions",
"when lambda is an lppm object",
"has changed (in spatstat 1.41-0 and later)",
"See help(linearKinhom)")
} else danger <- TRUE
lambdaX <-
if(is.vector(lambda)) lambda else
if(is.function(lambda)) lambda(XX$x, XX$y, ...) else
if(is.im(lambda)) safelookup(lambda, XX) else
if(inherits(lambda, "linim")) safelookup(as.im(lambda), XX) else
if(is.ppm(lambda) || inherits(lambda, "lppm"))
predict(lambda, locations=as.data.frame(XX)) else
stop(paste(lambdaname, "should be",
"a numeric vector, function, pixel image, or fitted model"))
if(!is.numeric(lambdaX))
stop(paste("Values of", lambdaname, "are not numeric"))
if((nv <- length(lambdaX)) != (np <- npoints(X)))
stop(paste("Obtained", nv, "values of", lambdaname,
"but point pattern contains", np, "points"))
if(any(lambdaX < 0))
stop(paste("Negative values of", lambdaname, "obtained"))
if(any(lambdaX == 0))
stop(paste("Zero values of", lambdaname, "obtained"))
if(danger)
attr(lambdaX, "dangerous") <- lambdaname
return(lambdaX)
}
linearKengine <- function(X, ..., r=NULL, reweight=NULL, denom=1,
correction="Ang", showworking=FALSE) {
# ensure distance information is present
X <- as.lpp(X, sparse=FALSE)
# extract info about pattern
np <- npoints(X)
# extract linear network
L <- domain(X)
W <- Window(L)
# determine r values
rmaxdefault <- 0.98 * boundingradius(L)
breaks <- handle.r.b.args(r, NULL, W, rmaxdefault=rmaxdefault)
r <- breaks$r
rmax <- breaks$max
#
type <- if(correction == "Ang") "L" else "net"
fname <- c("K", type)
ylab <- substitute(K[type](r), list(type=type))
#
if(np < 2) {
# no pairs to count: return zero function
zeroes <- numeric(length(r))
df <- data.frame(r = r, est = zeroes)
K <- fv(df, "r", ylab,
"est", . ~ r, c(0, rmax),
c("r", makefvlabel(NULL, "hat", fname)),
c("distance argument r", "estimated %s"),
fname = fname)
return(K)
}
# compute pairwise distances
D <- pairdist(X)
#--- compile into K function ---
if(correction == "none" && is.null(reweight)) {
# no weights (Okabe-Yamada)
K <- compileK(D, r, denom=denom, fname=fname)
K <- rebadge.fv(K, ylab, fname)
unitname(K) <- unitname(X)
return(K)
}
if(correction == "none")
edgewt <- 1
else {
# inverse m weights (Wei's correction)
# determine tolerance
toler <- default.linnet.tolerance(L)
# compute m[i,j]
m <- matrix(1, np, np)
for(j in 1:np)
m[ -j, j] <- countends(L, X[-j], D[-j,j], toler=toler)
if(any(uhoh <- (m == 0))) {
warning("Internal error: disc boundary count equal to zero")
m[uhoh] <- 1
}
edgewt <- 1/m
}
# compute K
wt <- if(!is.null(reweight)) edgewt * reweight else edgewt
K <- compileK(D, r, weights=wt, denom=denom, fname=fname)
# tack on theoretical value
K <- bind.fv(K, data.frame(theo=r),
makefvlabel(NULL, NULL, fname, "theo"),
"theoretical Poisson %s")
K <- rebadge.fv(K, ylab, fname)
unitname(K) <- unitname(X)
fvnames(K, ".") <- rev(fvnames(K, "."))
# show working
if(showworking)
attr(K, "working") <- list(D=D, wt=wt)
attr(K, "correction") <- correction
return(K)
}