Revision 9b559ffb328c82bf9b07288f5e1ba1ee9610e23c authored by Adrian Baddeley on 23 August 2009, 15:40:20 UTC, committed by cran-robot on 23 August 2009, 15:40:20 UTC
1 parent 4b1b757
Kmulti.S
#
# Kmulti.S
#
# Compute estimates of cross-type K functions
# for multitype point patterns
#
# $Revision: 5.30 $ $Date: 2009/07/24 06:08:22 $
#
#
# -------- functions ----------------------------------------
# Kcross() cross-type K function K_{ij}
# between types i and j
#
# Kdot() K_{i\bullet}
# between type i and all points regardless of type
#
# Kmulti() (generic)
#
#
# -------- standard arguments ------------------------------
# X point pattern (of class 'ppp')
# including 'marks' vector
# r distance values at which to compute K
#
# -------- standard output ------------------------------
# A data frame with columns named
#
# r: same as input
#
# trans: K function estimated by translation correction
#
# iso: K function estimated by Ripley isotropic correction
#
# theo: K function for Poisson ( = pi * r ^2 )
#
# border: K function estimated by border method
# using standard formula (denominator = count of points)
#
# bord.modif: K function estimated by border method
# using modified formula
# (denominator = area of eroded window
#
# ------------------------------------------------------------------------
"Lcross" <- function(X, i, j, ...) {
if(missing(i)) i <- levels(marks(X))[1]
if(missing(j)) j <- levels(marks(X))[2]
K <- Kcross(X, i, j, ...)
L <- eval.fv(sqrt(K/pi))
# relabel the fv object
L <- rebadge.fv(L,
substitute(L[i,j](r),
list(i=paste(i),j=paste(j))),
"Lcross",
new.yexp=substitute(L[list(i,j)](r),
list(i=paste(i),j=paste(j))))
return(L)
}
"Ldot" <- function(X, i, ...) {
if(missing(i)) i <- levels(marks(X))[1]
K <- Kdot(X, i, ...)
L <- eval.fv(sqrt(K/pi))
# relabel the fv object
L <- rebadge.fv(L,
substitute(Ldot[i](r), list(i=paste(i))),
"Ldot")
return(L)
}
"Kcross" <-
function(X, i, j, r=NULL, breaks=NULL,
correction =c("border", "isotropic", "Ripley", "translate") , ...)
{
verifyclass(X, "ppp")
if(!is.marked(X))
stop("point pattern has no marks (no component 'marks')")
if(missing(correction))
correction <- NULL
marx <- marks(X)
if(missing(i))
i <- levels(marx)[1]
if(missing(j))
j <- levels(marx)[2]
I <- (marx == i)
if(!any(I))
stop(paste("No points have mark i =", i))
if(i == j) {
result <- Kest(X[I],
r=r, breaks=breaks, correction=correction, ...)
} else {
J <- (marx == j)
if(!any(J))
stop(paste("No points have mark j =", j))
result <- Kmulti(X, I, J,
r=r, breaks=breaks, correction=correction, ...)
}
result <-
rebadge.fv(result,
substitute(Kcross[i,j](r), list(i=paste(i),j=paste(j))),
"Kcross",
new.yexp=substitute(K[list(i,j)](r),
list(i=paste(i),j=paste(j))))
return(result)
}
"Kdot" <-
function(X, i, r=NULL, breaks=NULL,
correction = c("border", "isotropic", "Ripley", "translate") , ...)
{
verifyclass(X, "ppp")
if(!is.marked(X))
stop("point pattern has no marks (no component 'marks')")
if(missing(correction))
correction <- NULL
marx <- marks(X)
if(missing(i))
i <- levels(marx)[1]
I <- (marx == i)
J <- rep(TRUE, X$n) # i.e. all points
if(!any(I)) stop(paste("No points have mark i =", i))
result <- Kmulti(X, I, J,
r=r, breaks=breaks, correction=correction, ...)
result <-
rebadge.fv(result,
substitute(Kdot[i](r), list(i=paste(i))),
"Kdot")
return(result)
}
"Kmulti"<-
function(X, I, J, r=NULL, breaks=NULL,
correction = c("border", "isotropic", "Ripley", "translate") , ...)
{
verifyclass(X, "ppp")
npoints <- X$n
x <- X$x
y <- X$y
W <- X$window
area <- area.owin(W)
correction.given <- !missing(correction) && !is.null(correction)
if(is.null(correction))
correction <- c("border", "isotropic", "Ripley", "translate")
correction <- pickoption("correction", correction,
c(none="none",
border="border",
"bord.modif"="bord.modif",
isotropic="isotropic",
Ripley="isotropic",
translate="translate",
best="best"),
multi=TRUE)
correction <- implemented.for.K(correction, W$type, correction.given)
if(!is.logical(I) || !is.logical(J))
stop("I and J must be logical vectors")
if(length(I) != npoints || length(J) != npoints)
stop(paste("The length of I and J must equal",
"the number of points in the pattern"))
if(!any(I)) stop("no points satisfy I")
if(!any(J)) stop("no points satisfy J")
nI <- sum(I)
nJ <- sum(J)
lambdaI <- nI/area
lambdaJ <- nJ/area
# r values
rmaxdefault <- rmax.rule("K", W, lambdaJ)
breaks <- handle.r.b.args(r, breaks, W, rmaxdefault=rmaxdefault)
r <- breaks$r
rmax <- breaks$max
# recommended range of r values
alim <- c(0, min(rmax, rmaxdefault))
# this will be the output data frame
# It will be given more columns later
K <- data.frame(r=r, theo= pi * r^2)
desc <- c("distance argument r", "theoretical Poisson %s")
K <- fv(K, "r", substitute(Kmulti(r), NULL),
"theo", , alim, c("r","%s[pois](r)"), desc, fname="Kmulti")
# find close pairs of points
XI <- X[I]
XJ <- X[J]
close <- crosspairs(XI, XJ, max(r))
# close$i and close$j are serial numbers in XI and XJ respectively;
# map them to original serial numbers in X
orig <- seq(npoints)
imap <- orig[I]
jmap <- orig[J]
iX <- imap[close$i]
jX <- jmap[close$j]
# eliminate any identical pairs
if(any(I & J)) {
ok <- (iX != jX)
if(!all(ok)) {
close$i <- close$i[ok]
close$j <- close$j[ok]
close$xi <- close$xi[ok]
close$yi <- close$yi[ok]
close$xj <- close$xj[ok]
close$yj <- close$yj[ok]
close$dx <- close$dx[ok]
close$dy <- close$dy[ok]
close$d <- close$d[ok]
}
}
# extract information for these pairs (relative to orderings of XI, XJ)
dcloseIJ <- close$d
icloseI <- close$i
jcloseJ <- close$j
# Compute estimates by each of the selected edge corrections.
if(any(correction == "none")) {
# uncorrected!
wh <- whist(dcloseIJ, breaks$val) # no weights
Kun <- cumsum(wh)/(lambdaI * lambdaJ * area)
K <- bind.fv(K, data.frame(un=Kun), "%s[un](r)",
"uncorrected estimate of %s",
"un")
}
if(any(correction == "border" | correction == "bord.modif")) {
# border method
# distance to boundary from each point of type I
bI <- bdist.points(XI)
# distance to boundary from first element of each (i, j) pair
bcloseI <- bI[icloseI]
# apply reduced sample algorithm
RS <- Kount(dcloseIJ, bcloseI, bI, breaks)
if(any(correction == "bord.modif")) {
denom.area <- eroded.areas(W, r)
Kbm <- RS$numerator/(denom.area * nI * nJ)
K <- bind.fv(K, data.frame(bord.modif=Kbm), "%s[bordm](r)",
"modified border-corrected estimate of %s",
"bord.modif")
}
if(any(correction == "border")) {
Kb <- RS$numerator/(lambdaJ * RS$denom.count)
K <- bind.fv(K, data.frame(border=Kb), "%s[bord](r)",
"border-corrected estimate of %s",
"border")
}
}
if(any(correction == "translate")) {
# translation correction
edgewt <- edge.Trans(XI[icloseI], XJ[jcloseJ], paired=TRUE)
wh <- whist(dcloseIJ, breaks$val, edgewt)
Ktrans <- cumsum(wh)/(lambdaI * lambdaJ * area)
rmax <- diameter(W)/2
Ktrans[r >= rmax] <- NA
K <- bind.fv(K, data.frame(trans=Ktrans), "%s[trans](r)",
"translation-corrected estimate of %s",
"trans")
}
if(any(correction == "isotropic")) {
# Ripley isotropic correction
edgewt <- edge.Ripley(XI[icloseI], matrix(dcloseIJ, ncol=1))
wh <- whist(dcloseIJ, breaks$val, edgewt)
Kiso <- cumsum(wh)/(lambdaI * lambdaJ * area)
rmax <- diameter(W)/2
Kiso[r >= rmax] <- NA
K <- bind.fv(K, data.frame(iso=Kiso), "%s[iso](r)",
"Ripley isotropic correction estimate of %s",
"iso")
}
# default is to display them all
attr(K, "fmla") <- . ~ r
unitname(K) <- unitname(X)
return(K)
}
Computing file changes ...