# R/linearpcf.R In spatstat.linnet: Linear Networks Functionality of the 'spatstat' Family

#### Documented in linearpcflinearpcfenginelinearpcfinhom

```#
# linearpcf.R
#
# \$Revision: 1.32 \$ \$Date: 2023/03/10 03:51:30 \$
#
# pair correlation function for point pattern on linear network
#
#
linearpcf <- function(X, r=NULL, ..., correction="Ang", ratio=FALSE) {
stopifnot(inherits(X, "lpp"))
correction <- pickoption("correction", correction,
c(none="none",
Ang="Ang",
best="Ang"),
multi=FALSE)
np <- npoints(X)
lengthL <- volume(domain(X))
# compute
samplesize <- npairs <- np * (np - 1)
denom <- npairs/lengthL
g <- linearpcfengine(X, r=r, ...,
denom=denom, correction=correction, ratio=ratio,
samplesize=samplesize)
# extract bandwidth
bw <- attr(g, "bw")
correction <- attr(g, "correction")
# set appropriate y axis label
switch(correction,
Ang  = {
ylab <- quote(g[L](r))
fname <- c("g", "L")
},
none = {
ylab <- quote(g[net](r))
fname <- c("g", "net")
})
# reattach bandwidth
attr(g, "bw") <- bw
attr(g, "correction") <- correction
return(g)
}

linearpcfinhom <- function(X, lambda=NULL, r=NULL,  ...,
correction="Ang", normalise=TRUE, normpower=1,
update=TRUE, leaveoneout=TRUE,
stopifnot(inherits(X, "lpp"))
loo.given <- !missing(leaveoneout)
correction <- pickoption("correction", correction,
c(none="none",
Ang="Ang",
best="Ang"),
multi=FALSE)
if(is.null(lambda))
warn.once("linearpcfinhomNULL",
"In linearpcfinhom the interpretation of 'lambda=NULL'",
"has changed (in spatstat.linnet 3.1 and later);",
"the function linearpcf is no longer invoked;",
"instead the intensity lambda is estimated by kernel smoothing")

if(normalise) {
check.1.real(normpower)
stopifnot(normpower >= 1)
}
lengthL <- volume(domain(X))
#
lambdaX <- resolve.lambda.lpp(X, lambda,
update=update,
leaveoneout=leaveoneout,
loo.given=loo.given,
sigma=sigma,
lambdaname="lambda",
#
invlam <- 1/lambdaX
invlam2 <- outer(invlam, invlam, "*")
denom <- if(!normalise) lengthL else
if(normpower == 1) sum(invlam) else
lengthL * (sum(invlam)/lengthL)^normpower
g <- linearpcfengine(X, ..., r=r,
reweight=invlam2, denom=denom,
correction=correction,
ratio=ratio)
# extract bandwidth
bw <- attr(g, "bw")
correction <- attr(g, "correction")
# set appropriate y axis label
switch(correction,
Ang  = {
ylab <- quote(g[L, inhom](r))
fname <- c("g", "list(L, inhom)")
},
none = {
ylab <- quote(g[net, inhom](r))
fname <- c("g", "list(net, inhom)")
})
# reattach bandwidth
attr(g, "bw") <- bw
attr(g, "correction") <- correction
attr(g, "dangerous") <- attr(lambdaX, "dangerous")
return(g)
}

linearpcfengine <- function(X, ..., r=NULL,
reweight=NULL,
denom=1, samplesize=NULL,
correction="Ang", ratio=FALSE) {
# ensure distance information is present
X <- as.lpp(X, sparse=FALSE)
np <- npoints(X)
# extract linear network
L <- domain(X)
W <- Window(L)
# determine r values
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("g", type)
ylab <- substitute(g[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)
g <- ratfv(df, NULL, 0,
"r", ylab,
"est", . ~ r, c(0, rmax),
c("r", makefvlabel(NULL, "hat", fname)),
c("distance argument r", "estimated %s"),
fname = fname,
ratio=ratio)
if(correction == "Ang") {
# tack on theoretical value
g <- bind.ratfv(g,
quotient = data.frame(theo=r),
denominator = 0,
labl = makefvlabel(NULL, NULL, fname, "theo"),
desc = "theoretical Poisson %s",
ratio=ratio)
}
attr(g, "correction") <- correction
return(g)
}
# compute pairwise distances
D <- pairdist(X)
#---  compile into pcf ---
if(correction == "none" && is.null(reweight)) {
g <- compilepcf(D, r, denom=denom, fname=fname, ratio=ratio,
samplesize=samplesize)
unitname(g) <- unitname(X)
attr(g, "correction") <- correction
return(g)
}
if(correction == "none") {
edgewt <- 1
} else {
## inverse m weights (Wei's correction)
## determine tolerance
toler <- default.linnet.tolerance(L)
## compute m[i,j]
m <- DoCountEnds(X, D, toler)
edgewt <- 1/m
}
# compute pcf
wt <- if(!is.null(reweight)) edgewt * reweight else edgewt
g <- compilepcf(D, r, weights=wt, denom=denom, ..., fname=fname, ratio=ratio,
samplesize=samplesize)
# extract bandwidth
bw <- attr(g, "bw")
# tack on theoretical value
g <- bind.ratfv(g,
quotient = data.frame(theo=rep.int(1,length(r))),
denominator = samplesize %orifnull% denom,
labl = makefvlabel(NULL, NULL, fname, "pois"),
desc = "theoretical Poisson %s",
ratio = ratio)
# tweak
unitname(g) <- unitname(X)
fvnames(g, ".") <- rev(fvnames(g, "."))
# tack on bandwidth again
attr(g, "bw") <- bw
attr(g, "correction") <- correction
return(g)
}
```

## Try the spatstat.linnet package in your browser

Any scripts or data that you put into this service are public.

spatstat.linnet documentation built on Nov. 2, 2023, 6:10 p.m.