knitr::opts_chunk$set(echo = TRUE) library(trajr) library(microbenchmark)
Mean vector characterising changes in direction:
# Using vector arithmetic. This is roughly twice as fast as using calculus angularMeanVectorArith <- function(angles) mean(complex(modulus = 1, argument = angles)) # Using calculus angularMeanCalculus <- function(angles) { phi <- atan2(sum(sin(angles)), sum(cos(angles))) r <- sqrt(sum(cos(angles)) ^ 2 + sum(sin(angles)) ^ 2) / length(angles) complex(modulus = r, argument = phi) } angles <- runif(10000, min = -pi, max = pi) microbenchmark(angularMeanVectorArith(angles), angularMeanCalculus(angles))
Empirical testing of Emax with generated trajectories.
set.seed(1) n <- 5 angularErrorSd <- runif(n, 0, 2) # Assess Emax for some different types of trajectories, all with the same angular error sd emaxr <- sapply(1:n, function(i) { TrajEmax(TrajGenerate(500, random = TRUE, angularErrorSd = angularErrorSd[i])) }) emaxd <- sapply(1:n, function(i) { TrajEmax(TrajGenerate(500, random = FALSE, angularErrorSd = angularErrorSd[i]), compass.direction = 0) }) emaxSm <- sapply(1:n, function(i) { TrajEmax(TrajRediscretize(TrajGenerate(500, random = TRUE, angularErrorSd = angularErrorSd[i]), R = 1)) }) emaxBig <- sapply(1:n, function(i) { TrajEmax(TrajRediscretize(TrajGenerate(500, random = TRUE, angularErrorSd = angularErrorSd[i]), R = 10)) }) emaxBig[emaxBig < 0] <- NA # can't take log of negative numbers par(mar = c(4, 4, 1, 1)) plot(rep(angularErrorSd, 4), c(emaxr, emaxd, emaxSm, emaxBig), log = 'xy', pch = '.', col = c(rep('red', n), rep('blue', n), rep('green', n), rep('black', n)), xlab = expression(sigma["angular error"]), ylab = expression("E"["max"])) labels <- c("Random walk", "Directed walk", "Random walk, step length 1", "Random walk, step length 10") legend("bottomleft", labels, pch = 16, col = c("red", "blue", "green", "black"), inset = 0.01)
Empirical testing of Sinuosity with generated trajectories.
set.seed(1) n <- 1000 angularErrorSd <- runif(n, 0, 2) # Assess sinuosity for some different types of trajectories, all with the same angular error sd sinur <- sapply(1:n, function(i) { TrajSinuosity(TrajGenerate(500, random = TRUE, angularErrorSd = angularErrorSd[i])) }) sinud <- sapply(1:n, function(i) { TrajSinuosity(TrajGenerate(500, random = FALSE, angularErrorSd = angularErrorSd[i]), compass.direction = 0) }) sinuSm <- sapply(1:n, function(i) { TrajSinuosity(TrajRediscretize(TrajGenerate(500, random = TRUE, angularErrorSd = angularErrorSd[i]), R = 1)) }) sinuBig <- sapply(1:n, function(i) { TrajSinuosity(TrajRediscretize(TrajGenerate(500, random = TRUE, angularErrorSd = angularErrorSd[i]), R = 10)) }) par(mar = c(4, 4, 1, 1)) plot(rep(angularErrorSd, 4), c(sinur, sinud, sinuSm, sinuBig), pch = '.', col = c(rep('red', n), rep('blue', n), rep('green', n), rep('black', n)), xlab = expression(sigma["angular error"]), ylab = "Sinuosity") legend("bottomright", c("Random walk", "Directed walk", "Random walk, step length 1", "Random walk, step length 10"), pch = 16, col = c("red", "blue", "green", "black"), inset = 0.01)
set.seed(1) n <- 1000 angularErrorSd <- runif(n, 0, pi) # Assess sinuosities for some different types of trajectories, all with the same angular error sd sins <- sapply(1:n, function(i) { c(TrajSinuosity(TrajGenerate(500, angularErrorSd = angularErrorSd[i])), TrajSinuosity2(TrajGenerate(500, angularErrorSd = angularErrorSd[i])) ) }) plot(sins[2,] ~ sins[1,], pch = 16, cex = .2, xlab = 'TrajSinuosity', ylab = 'TrajSinuosity2') abline(a = 0, b = 1, col = "red")
set.seed(1) n <- 1000 angularErrorSd <- runif(n, 0, 2) xy <- sapply(1:n, function(i) { trj <- TrajGenerate(500, angularErrorSd = angularErrorSd[i]) c(Emax = TrajEmax(trj), sinuosity = TrajSinuosity(trj)) }) xy <- xy[, which(xy[1,] > 0)] plot(log(xy[1,]), xy[2,], pch = '.', xlab = expression(log(E[max])), ylab = "Sinuosity")
set.seed(1) n <- 20 angularErrorSd <- runif(n, 0, 2) stepCount <- round(runif(n, 20, 100)) # Use the same step sizes for all trajectories #trj <- TrajGenerate((max(stepCount) + min(stepCount)) / 2) stepSizes <- exp(log(10) * seq(log10(.1), log10(7), length.out = 20)) dim <- sapply(1:n, function(i) { trj <- TrajGenerate(stepCount[i], random = TRUE, angularErrorSd = angularErrorSd[i]) c(TrajFractalDimension(trj, stepSizes), TrajFractalDimension(trj, stepSizes, FALSE, FALSE)) }) plot(angularErrorSd, dim[1,], pch = 16, col = rgb((stepCount - min(stepCount)) / (max(stepCount) - min(stepCount)), 0.2, 0.1), ylim = range(dim), xlab = expression(sigma["angular error"]), ylab = "Fractal dimension") col <- rgb(0.2, (stepCount - min(stepCount)) / (max(stepCount) - min(stepCount)), .01) points(angularErrorSd, dim[2,], pch = 3, col = col)
Plot of path length to divider size for a range of divider sizes, Levy flight and correlated random walk.
set.seed(1) .plotFractalValues <- function(trj, ...) { muL <- mean(TrajStepLengths(trj)) stepSizes <- TrajLogSequence(0.5 * muL, 5 * muL, 20) v <- TrajFractalDimensionValues(trj, stepSizes) d <- TrajFractalDimension(trj, stepSizes) l <- stats::lm(pathlength ~ stepsize, data = v) plot(v, log = "xy", pch = 16, ...) #abline(l, untf = TRUE) mtext(sprintf("Fractal dimension = %g", d), 1, -1.1, adj = 0.05) } lf <- TrajGenerate(linearErrorDist = rcauchy) crw <- TrajGenerate() par(mfrow = c(2, 1)) .plotFractalValues(lf, main = "Levy flight") .plotFractalValues(crw, main = "Correlated random walk", col = "red")
Following is an implementation of a simulation study from (Cheung at el., 2007). It plots the expected endpoints of multiple random trajectories after different numbers of steps, illustrating the difference in expected distance travelled between directed and random walks. It demonstrates that, for travel without a compass (i.e. an external means to determine direction), the distance travelled along an intended direction does not increase in proportion to the number of steps, which implies an upper limit to the distance that can be travelled without a compass.
set.seed(1) par(mfrow = c(2, 3), mar = c(4, 4, 3, 1)) ntraj <- 1000 # Cheung at el. used 20000, but it's slow # Define a function to simulate journeys of different lengths simulateJourneys <- function(random, nTrajectories = ntraj) { generate <- function(i) { # These parameters replicate the simulations in (Cheung et al., 2007) TrajGenerate(n = 500, random = random, stepLength = 2, linearErrorSd = 0.2, angularErrorSd = 0.5) } trjs <- lapply(1:nTrajectories, generate) # Use same plot axes as in the original article if (random) { xlims <- list(c(-20, 20), c(-20, 20), c(-20, 20), c(-40, 40), c(-170, 170), c(-400, 400)) ylims <- xlims } else { xlims <- list(c(-20, 20), c(-20, 20), c(-20, 20), c(20, 50), c(140, 200), c(810, 950)) ylims <- list(c(-20, 20), c(-20, 20), c(-20, 20), c(-20, 20), c(-30, 30), c(-70, 70)) } # 6 plots in 2 rows par(mfrow = c(2, 3)) steps <- c(1, 2, 5, 20, 100, 500) for(i in 1:length(steps)) { step <- steps[i] title <- sprintf("%d %s", step, ifelse(step == 1, "step", "steps")) # Plot the end-points of each journey after "step" steps endPoints <- t(sapply(trjs, function(t) t[step + 1, c('x', 'y')])) plot(endPoints, pch = '.', asp = 1, main = title, xlim = xlims[[i]], ylim = ylims[[i]]) } } # Call the function simulateJourneys(random = FALSE)
set.seed(1) par(mfrow = c(2, 3), mar = c(4, 4, 3, 1)) simulateJourneys(random = TRUE)
Benhamou (2006) defines a method, called backward pattern analysis, to determine whether a trajectory is a random or directed walk. Here we test it on a number of randomly generated trajectories, and report the failure rate.
# Returns TRUE if trj is a random trajectory according to backward pattern analysis (Benhamou, 2006) isRandom <- function(trj) { # For each pair of points along the trajectory, # calculate distance between the points D <- sapply(1:(nrow(trj) - 1), function(j) TrajDistance(trj, j, nrow(trj))) # calculate path length between points L <- sapply(1:(nrow(trj) - 1), function(j) TrajLength(trj, j, nrow(trj))) # Calculate linear regression with 0 intercept: D vs L dir <- lm(D ~ L - 1) # D vs sqrt(L) sL <- sqrt(L) rnd <- lm(D ~ sL - 1) # Pick the linear regression with the smallest sum of squared residuals sum(rnd$residuals ^ 2) < sum(dir$residuals ^ 2) } # Generate some random trajectories, randomly choosing # between random search path and directed walks set.seed(1) search <- sample(c(TRUE,FALSE), 1000, TRUE) # search/directed flag results <- sapply(1:length(search), function(i) { c(search[i], isRandom(TrajGenerate(100, random = search[i]))) }) # Flags whether the algorithm correctly identified the trajectory type correct <- results[1,] == results[2,] t <- table(data.frame(correct, results[1,])) cat(sprintf("Backward pattern analysis misclassified %g%% of random walks, and %g%% of directed walks\n", 100 * t[1,2] / sum(t), 100 * t[1,1] / sum(t)))
Following is a plot of the values used to calculate the fractal dimension of an approximation of an actual fractal curve. The points should lie on a straight line if the curve is a fractal. Also shows an inflexion point (i.e. horizontal line, then sloped line) since the value of n used to create the curve defines the scale down to which it is fractal.
# Adapted from http://www.soc.napier.ac.uk/~andrew/hilbert.html Hilbert <- function(x0, y0, xis, xjs, yis, yjs, n) { # Preallocate a matrix result <- matrix(NA, 4 ^ n, 2) row <- 1 .h <- function(x0, y0, xis, xjs, yis, yjs, n) { # x0 and y0 are the coordinates of the bottom left corner # xis & xjs are the i & j components of the unit x vector this frame # similarly yis and yjs if (n <= 0) { result[row, 1] <<- x0+(xis+yis)/2 result[row, 2] <<- y0+(xjs+yjs)/2 row <<- row + 1 } else { .h(x0, y0, yis/2, yjs/2, xis/2, xjs/2, n-1) .h(x0+xis/2, y0+xjs/2 ,xis/2, xjs/2, yis/2, yjs/2, n-1) .h(x0+xis/2+yis/2, y0+xjs/2+yjs/2, xis/2, xjs/2, yis/2, yjs/2,n-1) .h(x0+xis/2+yis, y0+xjs/2+yjs, -yis/2,-yjs/2, -xis/2, -xjs/2,n-1) } } .h(x0, y0, xis, xjs, yis, yjs, n) result } trj <- TrajFromCoords(as.data.frame(Hilbert(0, 0, 300, 0, 0, 300, 8))) plot(TrajFractalDimensionValues(trj, TrajLogSequence(.1, 400, 20)), log = 'xy', pch = 16) trj <- TrajFromCoords(as.data.frame(Hilbert(0, 0, 300, 0, 0, 300, 4))) points(TrajFractalDimensionValues(trj, TrajLogSequence(.1, 300, 40)), pch = 16, col = "red")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.