random.S
#
# random.S
#
# Functions for generating random point patterns
#
# $Revision: 4.23 $ $Date: 2006/10/17 09:08:09 $
#
#
# runifpoint() n i.i.d. uniform random points ("binomial process")
#
# runifpoispp() uniform Poisson point process
#
# rpoispp() general Poisson point process (thinning method)
#
# rpoint() n independent random points (rejection/pixel list)
#
# rMaternI() Mat'ern model I
# rMaternII() Mat'ern model II
# rSSI() Simple Sequential Inhibition process
#
# rNeymanScott() Neyman-Scott process (generic)
# rMatClust() Mat'ern cluster process
# rThomas() Thomas process
#
# rthin() independent random thinning
#
# Examples:
# u01 <- owin(0:1,0:1)
# plot(runifpoispp(100, u01))
# X <- rpoispp(function(x,y) {100 * (1-x/2)}, 100, u01)
# X <- rpoispp(function(x,y) {ifelse(x < 0.5, 100, 20)}, 100)
# plot(X)
# plot(rMaternI(100, 0.02))
# plot(rMaternII(100, 0.05))
#
"runifrect" <-
function(n, win=owin(c(0,1),c(0,1)))
{
# no checking
x <- runif(n, min=win$xrange[1], max=win$xrange[2])
y <- runif(n, min=win$yrange[1], max=win$yrange[2])
return(ppp(x, y, window=win))
}
"runifdisc" <-
function(n, radius=1, centre=c(0,0), ...)
{
# i.i.d. uniform points in the disc of radius r and centre (x,y)
disque <- disc(centre=centre, radius=radius, ...)
theta <- runif(n, min=0, max= 2 * pi)
s <- sqrt(runif(n, min=0, max=radius^2))
return(ppp(centre[1] + s * cos(theta),
centre[2] + s * sin(theta),
window=disque))
}
"runifpoint" <-
function(n, win=owin(c(0,1),c(0,1)), giveup=1000)
{
win <- as.owin(win)
if(n == 0)
return(ppp(numeric(0), numeric(0), window=win))
switch(win$type,
rectangle = {
return(runifrect(n, win))
},
mask = {
dx <- win$xstep
dy <- win$ystep
# extract pixel coordinates and probabilities
xpix <- as.vector(raster.x(win)[win$m])
ypix <- as.vector(raster.y(win)[win$m])
# select pixels with equal probability
id <- sample(seq(xpix), n, replace=TRUE)
# extract pixel centres and randomise within pixels
x <- xpix[id] + runif(n, min= -dx/2, max=dx/2)
y <- ypix[id] + runif(n, min= -dy/2, max=dy/2)
return(ppp(x, y, window=win))
},
polygonal={
# rejection method
# initialise empty pattern
x <- numeric(0)
y <- numeric(0)
X <- ppp(x, y, window=win)
#
# rectangle in which trial points will be generated
box <- bounding.box(win)
#
ntries <- 0
repeat {
ntries <- ntries + 1
# generate trial points in batches of n
qq <- runifrect(n, box)
# retain those which are inside 'win'
qq <- qq[win]
# add them to result
X <- superimpose(X, qq)
# if we have enough points, exit
if(X$n > n)
return(X[1:n])
else if(X$n == n)
return(X)
# otherwise get bored eventually
else if(ntries >= giveup)
stop(paste("Gave up after", giveup * n, "trials,",
np, "points accepted"))
}
})
stop("Unrecognised window type")
}
"runifpoispp" <-
function(lambda, win = owin(c(0,1),c(0,1))) {
win <- as.owin(win)
if(!is.numeric(lambda) || length(lambda) > 1 || lambda < 0)
stop("Intensity lambda must be a single number >= 0")
if(lambda == 0) # return empty pattern
return(ppp(numeric(0), numeric(0), window=win))
# generate Poisson process in enclosing rectangle
box <- bounding.box(win)
mean <- lambda * area.owin(box)
n <- rpois(1, mean)
X <- runifpoint(n, box)
# trim to window
if(win$type != "rectangle")
X <- X[win]
return(X)
}
rpoint <- function(n, f, fmax=NULL,
win=unit.square(), ..., giveup=1000,verbose=FALSE) {
if(missing(f) || (is.numeric(f) && length(f) == 1))
# uniform distribution
return(runifpoint(n, win, giveup))
# non-uniform distribution....
if(!is.function(f) && !is.im(f))
stop(paste(sQuote("f"),
"must be either a function or an",
sQuote("im"), "object"))
if(is.im(f)) {
# ------------ PIXEL IMAGE ---------------------
wf <- as.owin(f)
if(n == 0)
return(ppp(numeric(0), numeric(0), window=wf))
w <- as.mask(wf)
M <- w$m
dx <- w$xstep
dy <- w$ystep
# extract pixel coordinates and probabilities
xpix <- as.vector(raster.x(w)[M])
ypix <- as.vector(raster.y(w)[M])
ppix <- as.vector(f$v[M]) # not normalised - OK
# select pixels
id <- sample(length(xpix), n, replace=TRUE, prob=ppix)
# extract pixel centres and randomise within pixels
x <- xpix[id] + runif(n, min= -dx/2, max=dx/2)
y <- ypix[id] + runif(n, min= -dy/2, max=dy/2)
return(ppp(x, y, window=wf))
}
# ------------ FUNCTION ---------------------
# Establish parameters for rejection method
verifyclass(win, "owin")
if(n == 0)
return(ppp(numeric(0), numeric(0), window=win))
if(is.null(fmax)) {
# compute approx maximum value of f
imag <- as.im(f, win, ...)
summ <- summary(imag)
fmax <- summ$max + 0.05 * diff(summ$range)
}
irregular <- (win$type != "rectangle")
box <- bounding.box(win)
X <- ppp(numeric(0), numeric(0), window=win)
ntries <- 0
# generate uniform random points in batches
# and apply the rejection method.
# Collect any points that are retained in X
repeat{
ntries <- ntries + 1
# proposal points
prop <- runifrect(n, box)
if(irregular)
prop <- prop[win]
if(prop$n > 0) {
fvalues <- f(prop$x, prop$y, ...)
paccept <- fvalues/fmax
u <- runif(prop$n)
# accepted points
Y <- prop[u < paccept]
if(Y$n > 0) {
# add to X
X <- superimpose(X, Y)
if(X$n >= n) {
# we have enough!
if(verbose)
cat(paste("acceptance rate = ",
round(100 * X$n/(ntries * n), 2), "%\n"))
return(X[1:n])
}
}
}
if(ntries > giveup)
stop(paste("Gave up after",giveup * n,"trials with",
X$n, "points accepted"))
}
invisible(NULL)
}
"rpoispp" <-
function(lambda, lmax=NULL, win = owin(c(0,1),c(0,1)), ...) {
# arguments:
# lambda intensity: constant, function(x,y,...) or image
# lmax maximum possible value of lambda(x,y,...)
# win default observation window (of class 'owin')
# ... arguments passed to lambda(x, y, ...)
win <- if(is.im(lambda))
rescue.rectangle(as.owin(lambda))
else
as.owin(win)
if(is.numeric(lambda))
# uniform Poisson
return(runifpoispp(lambda, win))
# inhomogeneous Poisson
# perform thinning of uniform Poisson
if(is.null(lmax)) {
imag <- as.im(lambda, win, ...)
summ <- summary(imag)
lmax <- summ$max + 0.05 * diff(summ$range)
}
if(is.function(lambda)) {
X <- runifpoispp(lmax, win) # includes sanity checks on `lmax'
if(X$n == 0) return(X)
prob <- lambda(X$x, X$y, ...)/lmax
u <- runif(X$n)
retain <- (u <= prob)
X <- X[retain, ]
return(X)
}
if(is.im(lambda)) {
X <- runifpoispp(lmax, win)
if(X$n == 0) return(X)
prob <- lambda[X]/lmax
u <- runif(X$n)
retain <- (u <= prob)
X <- X[retain, ]
return(X)
}
stop(paste(sQuote("lambda"), "must be a constant, a function or an image"))
}
"rMaternI" <-
function(kappa, r, win = owin(c(0,1),c(0,1)))
{
X <- rpoispp(kappa, win=win)
if(X$n <= 1) return(X)
d <- nndist(X$x, X$y)
qq <- X[d > r]
return(qq)
}
"rMaternII" <-
function(kappa, r, win = owin(c(0,1),c(0,1)))
{
X <- rpoispp(kappa, win=win)
if(X$n <= 1) return(X)
# matrix of pairwise distances
d <- pairdist(X$x, X$y)
close <- (d <= r)
# random order 1:n
age <- sample(seq(X$n), X$n, replace=FALSE)
earlier <- outer(age, age, ">")
conflict <- close & earlier
# delete <- apply(conflict, 1, any)
delete <- matrowany(conflict)
qq <- X[ !delete]
return(qq)
}
"rSSI" <-
function(r, n, win = owin(c(0,1),c(0,1)), giveup = 1000)
{
# Simple Sequential Inhibition process
# fixed number of points
# Naive implementation, proposals are uniform
win <- as.owin(win)
X <- ppp(numeric(0),numeric(0), window=win)
r2 <- r^2
if(n * pi * r2/4 > area.owin(win))
stop(paste("Window is too small to fit", n, "points",
"at minimum separation", r))
ntries <- 0
while(ntries < giveup) {
ntries <- ntries + 1
qq <- runifpoint(1, win)
x <- qq$x[1]
y <- qq$y[1]
if(X$n == 0 || all(((x - X$x)^2 + (y - X$y)^2) > r2))
X <- superimpose(X, qq)
if(X$n == n)
return(X)
}
warning(paste("Gave up after", giveup,
"attempts with only", X$n, "points placed out of", n))
return(X)
}
"rNeymanScott" <-
function(kappa, rmax, rcluster, win = owin(c(0,1),c(0,1)), ..., lmax=NULL)
{
# Generic Neyman-Scott process
# Implementation for bounded cluster radius
#
# 'rcluster' is a function(x,y) that takes the coordinates
# (x,y) of the parent point and generates a list(x,y) of offspring
#
# "..." are arguments to be passed to 'rcluster()'
#
win <- as.owin(win)
# Generate parents in dilated window
frame <- bounding.box(win)
dilated <- owin(frame$xrange + c(-rmax, rmax),
frame$yrange + c(-rmax, rmax))
if(is.im(kappa) && !is.subset.owin(as.owin(kappa), dilated))
stop(paste("The window in which the image",
sQuote("kappa"),
"is defined\n",
"is not large enough to contain the dilation of the window",
sQuote("win")))
parents <- rpoispp(kappa, lmax=lmax, win=dilated)
#
# initialise offspring pattern and offspring-to-parent map
result <- ppp(numeric(0), numeric(0), window = win)
parentid <- numeric(0)
# generate clusters
if(parents$n > 0) {
for(i in seq(parents$n)) {
# generate random offspring of i-th parent point
cluster <- rcluster(parents$x[i], parents$y[i], ...)
cluster <- ppp(cluster$x, cluster$y, window=frame)
# trim to window
cluster <- cluster[win]
# add to pattern
result <- superimpose(result, cluster)
# update offspring-to-parent map
parentid <- c(parentid, rep(i, cluster$n))
}
}
attr(result, "parents") <- parents
attr(result, "parentid") <- parentid
return(result)
}
"rMatClust" <-
function(kappa, r, mu, win = owin(c(0,1),c(0,1)))
{
# Matern Cluster Process with Poisson (mu) offspring distribution
#
poisclus <- function(x0, y0, radius, mu) {
n <- rpois(1, mu)
return(runifdisc(n, centre=c(x0, y0),
radius=radius))
}
result <- rNeymanScott(kappa, r, poisclus,
win, radius=r, mu=mu)
return(result)
}
"rThomas" <-
function(kappa, sigma, mu, win = owin(c(0,1),c(0,1)))
{
# Thomas process with Poisson(mu) number of offspring
# at isotropic Normal(0,sigma^2) displacements from parent
#
thomclus <- function(x0, y0, sigma, mu) {
n <- rpois(1, mu)
x <- rnorm(n, mean=x0, sd=sigma)
y <- rnorm(n, mean=y0, sd=sigma)
return(list(x=x, y=y))
}
result <- rNeymanScott(kappa, 4 * sigma, thomclus,
win, sigma=sigma, mu=mu)
return(result)
}
rstrat <- function(win=square(1), nx, ny=nx, k=1) {
win <- as.owin(win)
stopifnot(nx >= 1 && ny >= 1)
stopifnot(k >= 1)
xy <- stratrand(win, nx, ny, k)
Xbox <- ppp(xy$x, xy$y, win$xrange, win$yrange)
X <- Xbox[win]
return(X)
}
rsyst <- function(win=square(1), nx, ny=nx, dx=NULL, dy=NULL) {
win <- as.owin(win)
xr <- win$xrange
yr <- win$yrange
if(!missing(nx)) {
stopifnot(nx >= 1 && ny >= 1)
if(!is.null(dx) || !is.null(dy))
stop("Do not give both nx and dx")
dx <- diff(xr)/nx
dy <- diff(yr)/ny
x0 <- seq(xr[1], xr[2], length=nx+1)
y0 <- seq(yr[1], yr[2], length=ny+1)
} else if(!is.null(dx)) {
stopifnot(dx > 0 & dy > 0)
if(!is.null(nx) || !is.null(ny))
stop("Do not give both nx and dx")
x0 <- seq(xr[1], xr[2], by=dx)
y0 <- seq(yr[1], yr[2], by=dy)
} else stop("Either (nx, ny) or (dx, dy) should be given")
xy0 <- expand.grid(x=x0, y=y0)
x <- xy0$x + runif(1, min = 0, max = dx)
y <- xy0$y + runif(1, min = 0, max = dy)
Xbox <- ppp(x, y, xr, yr)
X <- Xbox[win]
return(X)
}
rcell <- function(win=square(1), nx, ny=nx, dx=NULL, dy=NULL) {
win <- as.owin(win)
xr <- win$xrange
yr <- win$yrange
if(!missing(nx)) {
stopifnot(nx >= 1 && ny >= 1)
if(!is.null(dx) || !is.null(dy))
stop("Do not give both nx and dx")
dx <- diff(xr)/nx
dy <- diff(yr)/ny
x0 <- seq(xr[1], xr[2], length=nx+1)
y0 <- seq(yr[1], yr[2], length=ny+1)
} else if(!is.null(dx)) {
stopifnot(dx > 0 & dy > 0)
if(!is.null(nx) || !is.null(ny))
stop("Do not give both nx and dx")
x0 <- seq(xr[1], xr[2], by=dx)
y0 <- seq(yr[1], yr[2], by=dy)
} else stop("Either (nx, ny) or (dx, dy) should be given")
mx <- length(x0)
my <- length(y0)
x <- numeric(0)
y <- numeric(0)
rcellnumber <- function(n) {
u <- runif(n, min=0, max=1)
k <- ifelse(u < 1/10, 0, ifelse(u < 89/90, 1, 10))
return(k)
}
for(ix in seq(mx))
for(iy in seq(my)) {
nij <- rcellnumber(1)
x <- c(x, x0[ix] + runif(nij, min=0, max=dx))
y <- c(y, y0[iy] + runif(nij, min=0, max=dy))
}
Xbox <- ppp(x, y, xr, yr)
X <- Xbox[win]
return(X)
}
rthin <- function(X, P, ...) {
verifyclass(X, "ppp")
if(is.numeric(P)) {
# vector of retention probabilities
pX <- P
if(length(pX) != X$n) {
if(length(pX) == 1)
pX <- rep(pX, X$n)
else
stop("Length of vector P does not match number of points of X")
}
if(any(is.na(pX)))
stop("P contains NA's")
} else if(is.function(P)) {
# function - evaluate it at points of X
pX <- P(X$x, X$y, ...)
if(length(pX) != X$n)
stop("Function P returned a vector of incorrect length")
if(!is.numeric(pX))
stop("Function P returned non-numeric values")
if(any(is.na(pX)))
stop("Function P returned some NA values")
prange <- range(pX)
} else if(is.im(P)) {
# image - look it up
if(!(P$type %in% c("integer", "real")))
stop("Values of image P should be numeric")
pX <- P[X, drop=FALSE]
if(any(is.na(pX)))
stop("some points of X lie outside the domain of image P")
} else
stop("Unrecognised format for P")
if(min(pX) < 0) stop("some probabilities are negative")
if(max(pX) > 1) stop("some probabilities are greater than 1")
retain <- (runif(length(pX)) < pX)
return(X[retain])
}