Simulation of point patterns

We now look at simulating point patterns, which can be useful to test theoretical properties and to compare with real data to be analysed.

Simulating Homogeneous Poisson Processes

In this section, we will look further into the simulation of homogeneous poisson processes. We will simulate a process with intensity $\lambda=100$ points per unit square.

Start by defining the following polygon and setting the value of the intensity:

oct <- owin(poly = list(x = c(0, 1/3, 2/3, 1, 1, 2/3,
    1/3, 0), y = c(1/3, 0, 0, 1/3, 2/3, 1, 1, 2/3)))
lambda <- 100
plot(oct)

Because it is easy to simulate a homogeneous Poisson process on a (unit) square, we will begin by using rejection sampling to generate the homogeneous process on the octagonal shape, oct.

Recall that the number of points inside the window will follow a Poisson distribution with mean $\lambda A$ , where $A$ is the area of the observation window. To compute the area of the observation window, we use,

A <- area.owin(oct)

So, we draw the number of points to simulate from a Poisson with mean lambda*A,

npts <- rpois(1, lambda = A * lambda)

So we will simulate 78 points. In the next chunk of code, we perform the rejection sampling, storing the accepted points in vectors x and y. Make sure you understand how this works.

x <- c()
y <- c()
while (length(x) <= npts) {
    xcand <- runif(1)
    ycand <- runif(1)
    test <- inside.owin(x = xcand, y = ycand, w = oct)
    if (test) {
        x <- c(x, xcand)
        y <- c(y, ycand)
    }
}
patt1 <- ppp(x = x, y = y, window = oct)

plot(patt1)
axis(1)
axis(2)

Exercise 1.6.

Using rejection sampling, simulate and plot data generated from a homogeneous Poisson process with intensity 50 on a triangle with vertices $(-1,-1)$ , $(3,1)$ and $(2,2)$ .

Exercise 1.7.

Produce similar realisations of homogeneous Poisson process on oct and tri using the built-in spatstat functions.

Simulating Inhomogeneous Poisson Processes

In this section, you will simulate an inhomogeneous Poisson processes with intensity,

\lambda(x_{1},x_{2})=10+100x_{1}+200x_{2},

on the window $[0,1]\times[0,2]$ . In an inhomogeneous process, the intensity is not constant over space; we will again use rejection sampling to simulate the process.

We begin by defining the a function that will return the intensity. The function will take a single vector argument, x.

lambda <- function(x) {
    return(10 + 100 * x[1] + 200 * x[2])
}

It is informative to first plot $\lambda(x)$ on a fine grid, to see how it is varying. Define the coordinates of the points where $\lambda(x)$ will be evaluated:

xseq <- seq(0, 1, length.out = 50)
yseq <- seq(0, 2, length.out = 100)
grid <- expand.grid(xseq, yseq)

Then we evaluate our function at each of these points

z <- apply(grid, 1, lambda)

A nice tool for visualising function values on a grid is the image.plot function from the package fields

zmat <- matrix(z, 50, 100)
library(fields)
image.plot(xseq, yseq, zmat, xlab = "x", ylab = "y",
    main = "lambda(x)")

Alternative visualisations include persp(xseq,yseq,zmat) and contour(xseq,yseq,zmat).

In order to proceed with rejection sampling, we need to know the highest value of the intensity function in the observation window. Considering the form of $\lambda$ , it is obvious that the maximum will be the value $\lambda_{\text{max}}=\lambda((1,2))$ .

lambdamax <- lambda(c(1, 2))
lambdamax

[1] 510

For the rejection scheme, we propose from a completely spatially random process (i.e. a Uniform density) on $W=[0,1]\times[0,2]$ . In a similar manner to the homogeneous example above, we will need to integrate $\lambda$ over $W$ in order to simulate the number of realisations to draw.

Exercise 1.8.

Show that the integral, $V$ , of $\lambda(x)$ over $W$ is equal to 520.

Hence the number of points we will simulate, npts, is

npts <- rpois(1, lambda = 520)

The main difference between this rejection algorithm and the last is that we keep each candidate

(x_{\text{cand}},y_{\text{cand}})

with probability

\frac{\lambda(x_{\text{cand}},y_{\text{cand}})}{\lambda_{\text{max}}}.

x <- c()
y <- c()
while (length(x) <= npts) {
    xcand <- runif(1)
    ycand <- runif(1, 0, 2)
    prob <- lambda(c(xcand, ycand))/lambdamax
    if (prob > runif(1)) {
        x <- c(x, xcand)
        y <- c(y, ycand)
    }
}
patt3 <- ppp(x = x, y = y, window = owin(xrange = c(0,
    1), yrange = c(0, 2)))

plot(patt3)
axis(1)
axis(2)

Make sure that you understand all of the steps carried out here.

Exercise 1.9.

Using a spatially-varying $\lambda(x)$ and observation window, $W$ , of your choice, simuate a realisation of an inhomogenegous Poisson process on $W$ with intensity $\lambda$ using rejection sampling. Check your answer by simulating a second realisation using the rpoispp function. Remember that you need to be able to integrate $\lambda$ over $W$ .

Remember to save your workspace before finishing the session.