Originally written around 2014. Updated 10.12.2022.
The package Kdirectional
is a collection of tools for analysing the isotropy of 2- and 3-dimensional point patterns. This document shows how to use some of these tools by analysing an synthetic example pattern.
In this vignette we use the following packages:
library(spatstat) library(Kdirectional)
library(spatstat) library(Kdirectional)
We simulate Strauss process and transform it. Assume the transformation is
comp <- 0.7 D <- diag( c(comp, 1/comp)) angle <- pi/8 R <- rotationMatrix3(az=angle)[-3,-3] # internal, not rgl one C <- R%*%D
that is, we compress the x-axis by factor 0.7 and stretch the y-axis to keep volume 1. Then we simulate the original process in a window that will be square after the transformation
set.seed(1) x0 <- rStrauss(500, gamma = 0.05, R = 0.03, W = affine(square(), solve(C)))
and transform it to the observed pattern
x <- affine(x0, C)
Check things went ok (define a little helper for plotting):
plotp <- function(...) plot(..., cex=.3) par(mfrow=c(1,2), mar=c(1,1,1,1)) plotp(x0) plotp(x)
The pattern on the left is what we want to estimate.
Check for starters the distribution of nearest neighbours vectors.
nna <- nnangle(coords(x)) ang <- nna[,1]-pi par(mfrow=c(1,2)) angles.flower(ang, ci = T, k=24, bty = "n") # Compare to spatstat rose(ang, unit = "rad", nclass = 24+1, main = "")
The blue rings give approximate 95% pointwise confidence intervals. Hard to see anything definite, perhaps some peakyness at the 2:30'oclock axis. Nearest neighbours are not very informative, so let's analyse the full set of pairwise vectors.
As in the paper [1], the Fry-points tell us interesting things about the isotropy.
f <- fry_ellipsoids(x, nvec=1:20, nangles = 30, eps=0.1, r_adjust = .2)
r_adjust=.2
reduces the computations to pairs with length within 20% window width, eps
adds a bit of angular smoothing (in radians), and nangles
chooses the amount of directions to look at.
Plot the estimated ellipsoids:
plot(f, zoom=.3, used_points = FALSE)
Oval shape is visible, some ambiguity of rotation. Check the ellipticity with a contrast $axis_x-axis_y=0$:
ci <- confint(f, nsim=500) print(summary(ci)[1:5,], digits=2) plot(ci, ylim = c(-1,1)*.1)
The first contour is poorly estimated as there are some short range pairs, kind of like Poisson noise. Others seem to hint of compression (the p-values and confidence intervals are biased and conservative).
Let's compute the average rotation, where we scale each ellipsoid to have det = 1 by setting just_rotation=TRUE
.
contours <- f$ellipsoids[-1] mean_e <- mean_ellipsoids(contours, just_rotation = TRUE) summary(mean_e)
Rotation seems to be non-zero (or non $\pm\pi$). Ignore the error variance, its not correct.
c(true=angle, est=mean_e$rot_angle + pi)
Backrotate the pattern:
Rhat <- rotationMatrix3(az=pi)[-3,-3]%*%mean_e$rot xr <- affine(x, solve(Rhat))
Next we estimate the compression by minimizing the difference of K-functions along the main axes. We need to provide a grid to optimise over:
gamma_vec <- seq(0.5, 1, by=0.05) grid <- cbind(gamma_vec, 1/gamma_vec)
Then compute the anisotropy profile:
r <- seq(0, 0.2, length=50) eps <- pi/10 ani <- anisotropy_profile(xr, grid = grid, r=r, eps=eps, power=2) ani0 <- anisotropy_profile(x, grid = grid, r=r, eps=eps, power=2)
We set the sector-K angle with eps
and use integral of squared differences (power=2
) as the distance.
Plot the profiles, scaled to 0-1 for plotting (it's the default):
plot(ani, type="b", scale=TRUE) lines(ani0, type="b", col=2)
( s <- summary(ani) )
Estimated compression
c(true=comp, est=s$estimate[1])
Backtransform:
Dhat <- s$estimate xhat <- affine(xr, solve(Dhat))
Recap of the steps:
par(mfrow=c(2,2), mar=c(1,1,1,1)) plotp(x) plotp(xr) plotp(xhat) plotp(x0)
Finally, let's illustrate the estimated transformation. Use simulation and estimation to get an ellipsoid object:
Chat <- Rhat %*% Dhat e0 <- ellipsoid_OLS( rellipsoid(1000, c(1,1), R = C) ) e <- ellipsoid_OLS( rellipsoid(1000, c(1,1), R = Chat) )
plot(e0, add=FALSE) plot(e, col=3)
[1] T Rajala, A Särkkä, C Redenbach, M Sormani (2015): Estimating geometric anisotropy in spatial point patterns, Spatial Statistics
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.