marked point processes code

Jonathan Monteleone (j.monteleone@auckland.ac.nz)
Wed, 28 Jan 1998 07:23:59 +1300


The following code computes the K functions
for multitype point patterns
(i.e. where the marks are discrete).

Contents:
Kcross.S estimation of cross K function K_{ij}(t)
between two types i and j
Kest.S estimation of K function K_{ii}(t)
within one type i

########################################################################
#
# Kcross.S Estimation of cross K function
# for bivariate point pattern
#
# $Revision: 1.1 $ $Date: 1997/09/05 17:06:30 $
#
# *** Experimental ****
#
# So far, implemented for rectangular window only.
# Border correction only.
#
#
# -------- functions ----------------------------------------
# Kcross.b1() compute estimate of K by border correction,
# denominator = count of points
#
# Kcross.b2() compute estimate of K
# denominator = area of eroded window
#
# crossdist() compute matrix of distances between
# each pair of data points
#
#
# -------- standard arguments ------------------------------
# x1,y1 coordinates of input data points (type 1)
# x2,y2 coordinates of input data points (type 2)
# xrange,yrange dimensions of observation window
# r distance values at which to compute K
#
# -------- example -----------------------------------------
#
# x1 <- runif(20)
# y1 <- runif(20)
# x2 <- runif(30)
# y2 <- runif(30)
# plot(c(0,1),c(0,1),type="n")
# points(x1,y1,pch=2)
# points(x2,y2,pch=3)
# r <- seq(0,0.25,length=50)
# k <- Kest.b2(x1,y1,x2,y2,r=r)
# plot(r, k, type="l")
# lines(r, pi * r^2, lty=2)
#
# ------------------------------------------------------------------------

"crossdist"<-
function(x1, y1, x2, y2)
{
if(length(x1) != length(y1))
stop("lengths of x1 and y1 do not match")
if(length(x2) != length(y2))
stop("lengths of x2 and y2 do not match")
n1 <- length(x1)
n2 <- length(x2)
X1 <- matrix(rep(x1, n2), ncol = n2)
Y1 <- matrix(rep(y1, n2), ncol = n2)
X2 <- matrix(rep(x2, n1), ncol = n1)
Y2 <- matrix(rep(y2, n1), ncol = n1)
d <- sqrt((X1 - t(X2))^2 + (Y1 - t(Y2))^2)
return(d)
}
"Kcross.b1"<-
function(x1, y1, x2, y2, xrange=c(0,1), yrange=c(0,1), r)
{
if(length(x1) != length(y1))
stop("lengths of x1 and y1 do not match")
if(length(x2) != length(y2))
stop("lengths of x2 and y2 do not match")
n1 <- length(x1)
n2 <- length(x2)
width <- abs(diff(xrange))
height <- abs(diff(yrange))
area <- width * height
lambda1 <- n1/area
lambda2 <- n2/area
d <- crossdist(x1, y1, x2, y2)
b <- bdist(x1, y1, xrange, yrange)
k <- rep(0, length(r))
for(i in 1:length(r)) {
close <- (d <= r[i])
nclose <- apply(close, 1, sum)
bok <- (b > r[i])
numerator <- sum(nclose[bok])
denominator <- sum(bok)
k[i] <- numerator/denominator
}
k <- k/lambda2
return(k)
}
"Kcross.b2"<-
function(x1, y1, x2, y2, xrange=c(0,1), yrange=c(0,1), r)
{
if(length(x1) != length(y1))
stop("lengths of x1 and y1 do not match")
if(length(x2) != length(y2))
stop("lengths of x2 and y2 do not match")
n1 <- length(x1)
n2 <- length(x2)
width <- abs(diff(xrange))
height <- abs(diff(yrange))
area <- width * height
lambda1 <- n1/area
lambda2 <- n2/area
d <- crossdist(x1, y1, x2, y2)
b <- bdist(x1, y1, xrange, yrange)
k <- rep(0, length(r))
for(i in 1:length(r)) {
close <- (d <= r[i])
nclose <- apply(close, 1, sum)
bok <- (b > r[i])
numerator <- sum(nclose[bok])
denominator <- (width - 2 * r[i]) * (height - 2 * r[i])
k[i] <- numerator/denominator
}
k <- k/(n1 * n2)
return(k)
}

#######################################################################
#
# Kest.S Estimation of K function
#
# $Revision: 1.1 $ $Date: 1997/08/08 09:38:37 $
#
# *** Experimental ****
#
# So far, implemented for rectangular window only.
# Border correction only.
#
#
# -------- functions ----------------------------------------
# Kest.b1() compute estimate of K by border correction,
# denominator = count of points
# Kest.b2() compute estimate of K
# denominator = area of eroded window
#
# bdist() compute vector of distances from
# each data point to boundary of window
#
# pairdist() compute matrix of distances between
# each pair of data points
#
#
# -------- standard arguments ------------------------------
# x,y coordinates of input data points
# xrange,yrange dimensions of observation window
# r distance values at which to compute K
#

"bdist"<-
function(x, y, xrange=c(0,1), yrange=c(0,1))
{
xd <- pmin(x - min(xrange), max(xrange) - x)
yd <- pmin(y - min(yrange), max(yrange) - y)
d <- pmin(xd, yd)
return(d)
}
"pairdist"<-
function(x, y)
{
if(length(x) != length(y))
stop("lengths of x and y do not match")
n <- length(x)
xx <- matrix(rep(x, n), nrow = n)
yy <- matrix(rep(y, n), nrow = n)
d <- sqrt((xx - t(xx))^2 + (yy - t(yy))^2)
return(d)
}
"Kest.b1"<-
function(x, y, xrange=c(0,1), yrange=c(0,1), r)
{
if(length(x) != length(y))
stop("lengths of x and y do not match")
npoints <- length(x)
width <- abs(diff(xrange))
height <- abs(diff(yrange))
area <- width * height
lambda <- npoints/area
d <- pairdist(x, y)
b <- bdist(x, y, xrange, yrange)
k <- rep(0, length(r))
for(i in 1:length(r)) {
close <- (d <= r[i])
nclose <- apply(close, 1, sum) - 1
bok <- (b > r[i])
numerator <- sum(nclose[bok])
denominator <- sum(bok)
k[i] <- numerator/denominator
}
k <- k/lambda
return(k)
}
"Kest.b2"<-
function(x, y, xrange=c(0,1), yrange=c(0,1), r)
{
if(length(x) != length(y))
stop("lengths of x and y do not match")
npoints <- length(x)
width <- abs(diff(xrange))
height <- abs(diff(yrange))
area <- width * height
lambda <- npoints/area
lambda2 <- (npoints * (npoints - 1))/(area^2)
d <- pairdist(x, y)
b <- bdist(x, y, xrange, yrange)
k <- rep(0, length(r))
for(i in 1:length(r)) {
close <- (d <= r[i])
nclose <- apply(close, 1, sum) - 1
bok <- (b > r[i])
numerator <- sum(nclose[bok])
denominator <- (width - 2 * r[i]) * (height - 2 * r[i])
k[i] <- numerator/denominator
}
k <- k/lambda2
return(k)
}