Objectives

  • Develop an intuition of the Brownian Bridge model using random walks.
  • Practice R Coding

General Idea

The central concept here is that if we have two known locations of an organism, for example two points from radiotelemetry data, we can ‘interpolate’ a set of possible random paths they took between the points.

We could simulate such paths using a random walk model.

However, there is an issue: while we can specify a starting point for a random walk, there’s no obviously way to force our random walk to go through a specific ending point. What to do?

Here is a possible workflow:

  1. Simulate a random walk originating at location 1.
  2. Rotate and stretch the path of the random walk so that it ends at location 2.
  3. Repeat the process many times to generate a probability surface.

Preliminaries

To create and plot the simulations, we’ll need some supporting functions.

Plot Function

We’ll re-use our random-walk plotting function from the in-class random walk exercise.

plot_sim = function(
    dat_sim,
    add = FALSE,
    ylab = "y", 
    xlab = "x",
    col_line = 1,
    col_start = 1, col_end = 2,
    cex_start = 2, cex_end = 2,
    pch_start = 16, pch_end = 16,
    ...)
{
  
  # If this is the first simulation, create a new plot
  if (!add)
  {
    plot(
      dat_sim[, 1], dat_sim[, 2],
      xlab = xlab, ylab = ylab,
      type = "l", asp = 1,
      col = col_line,
      ...)
  }
  
  # If adding to an existing plot:
  if (add)
  {
    lines(
      dat_sim[, 1], dat_sim[, 2],
      col = col_line,
      ...)
  }
  
  points(
    dat_sim[1, 1], dat_sim[1, 2],
    pch = pch_start, col = col_start,
    cex = cex_start)
  points(
    dat_sim[nrow(dat_sim), 1],
    dat_sim[nrow(dat_sim), 2],
    pch = pch_end, col = col_end,
    cex = cex_end)
}

Correlated Random Walk Function

We can use my implementation of the correlated random walk function:

r_walk_correlated = function(
    n_steps,
    step_min = 0.1,
    step_max = 1,
    cone_width = pi/4,
    initial_direction = 0,
    x0 = 0, y0 = 0)
{
  # Empty matrix to keep track of the position
  data_steps = matrix(0, ncol = 2, nrow = n_steps + 1)
  
  # before step, set x1 y1 to x0, y0
  x_i = x0; y_i = y0
  
  direction_i = initial_direction
  
  # Loop for the walk:
  for(i in 1:n_steps)
  {
    dir_add = runif(1, min = -cone_width/2, max = cone_width/2)
    direction_i = direction_i + dir_add
    
    step_length = runif(n = 1, min = step_min, max = step_max)
    
    # horizontal
    x_i = x_i + cos(direction_i) * step_length
    # vertical
    y_i = y_i + sin(direction_i) * step_length
    
    data_steps[i+1, 1] = x_i
    data_steps[i+1, 2] = y_i
  }
  return(data_steps)
}

Helper Functions

We also need a few supporting functions to perform the following tasks:

  1. Calculate the angle (in radians) between two coordinates.
  2. Rotate a set of coordinates by a specified angle.

Here are some possible implementations:

Rotation Function

Rotation function using a rotation matrix

# Rotate a single coordinate
rot = function(coord, theta)
{
  mat = matrix(
    c(cos(theta), -sin(theta),
      sin(theta),  cos(theta)),
    nrow = 2, byrow = T
  )
  return(mat %*% coord)
}

# rotate a set of coordinates
rot_coords = function(coords, theta)
{
  t(apply(coords, 1, function(x) rot(matrix(x), theta)))
}

Angle Function

Angle function using the inverse tangent function:

get_angle = function(coord1, coord2)
{
  delta_x = coord1[1] - coord2[1]
  delta_y = coord1[2] - coord2[2]
  atan(delta_y / delta_x)
}

Scale Function

Scale a trajectory according to its distance in the x-direction:

stretch_x = function(coords, length_x = 1)
{
  x0 = head(coords, 1)[1]
  x1 = tail(coords, 1)[1]
  x_stretch = x1 - x0
  stretch_factor = length_x / x_stretch
  return(coords * stretch_factor)
  # return(coords / x_stretch)
}

Helper Function Tests

Let’s see our helper functions in action:

Simulation data

set.seed(11)
sim_coords = r_walk_correlated(10)
plot_sim(sim_coords)

Now rotate by -45, 90, and 270 degrees:

plot_sim(sim_coords, xlim = c(-6, 6))
plot_sim(rot_coords(sim_coords, pi / 4), xlim = c(), add = T)
plot_sim(rot_coords(sim_coords, pi / 2), xlim = c(), add = T)
plot_sim(rot_coords(sim_coords, (3 * pi) / 2), xlim = c(), add = T)

Calculate the angle between the starting point and end point:

get_angle(sim_coords[1,], tail(sim_coords, 1))
## [1] -0.5014947

Rotate path so ending coordinate is on x-axis:

sim_coords1 = rot_coords(
  sim_coords,
  -get_angle(sim_coords[1,], tail(sim_coords, 1))
)
plot_sim(sim_coords, xlim = c(-6, 6))
plot_sim(sim_coords1, xlim = c(-6, 6), add = T)
abline(h = 0, lty = 2)

Scale the path so that it goes between the origin and (1, 0)

sim_coords2 = stretch_x(sim_coords, length_x = 1)
sim_coords3 = stretch_x(sim_coords1, length_x = 1)
plot_sim(sim_coords, xlim = c(0, 6))
plot_sim(sim_coords1, add = T)
plot_sim(sim_coords2, add = T)
plot_sim(sim_coords3, add = T)
abline(h = 0, lty = 2)
abline(v = 1, lty = 2)

Everything Function

Now we need to combine everything in to a master function:

rotate_scale = function(coords, length_x = 1)
{
  coords_rotated = rot_coords(
    coords,
    -get_angle(coords[1,], tail(coords, 1))
  )
  coords_scaled = stretch_x(coords_rotated, length_x)
  return(coords_scaled)
}

plot_sim(rotate_scale(sim_coords))

Brownian Path Simulation

Now we’re ready to run some simulations!

First plot a path:

n_steps = 10
n_sims = 1000

{
  plot_sim(
    rotate_scale(r_walk_correlated(n_steps)),
    xlim = c(-0.2, 1.2))
  
  for (i in 1:n_sims)
    plot_sim(
      rotate_scale(r_walk_correlated(n_steps)),
      add = T,
      col_line = adjustcolor("steelblue", 0.1))
}

Success!

Now try out the simulation with different settings for the random walk function.
Some things to try:

  • Adjust the cone width
  • Try different numbers of steps
  • Try different numbers of simulations