R/trefoil.r

Defines functions jd_trefoil rs_trefoil sample_trefoil

Documented in sample_trefoil

#' @title Sample (with noise) from trefoil knot
#'
#' @description These functions generate uniform samples from trefoil knot in
#'   3-dimensional space, optionally with noise.
#'
#' @details The trefoil knot is the simplest nontrivial knot and contains three
#'   unique crossings in three dimensional space. This uniform sample is
#'   generated by rejection sampling. This process allows for simulation of
#'   random samples from the trefoil knot distribution, by using random samples
#'   from a more convenient distribution. It applies rejection/acceptance
#'   criterion such that the samples that are accepted are to be distributed as
#'   if they were from the target distribution. The uniform sample is generated
#'   through a rejection sampling process as described by Diaconis, Holmes, and
#'   Shahshahani (2013).

#' @template ref-diaconis2013
#' 

#' @name trefoil
#' @param n Number of observations.
#' @param sd Standard deviation of (independent multivariate) Gaussian noise.
#' @example inst/examples/ex-trefoil.r

NULL

#' @rdname trefoil
#' @export

sample_trefoil <- function(n, sd = 0){
  theta <- rs_trefoil(n)
  # Applies modified theta values to parametric equation of trefoil knot
  res <- cbind(
    x = sin(theta) + 2*sin(2*theta),
    y = cos(theta) - 2*cos(2*theta),
    z = -sin(3*theta)
  )
  #Adds Gaussian noise to figure 8
  add_noise(res, sd = sd)
}

# Rejection sampler
rs_trefoil <- function(n){
  x <- c()
  # Keep looping until desired number of observations is achieved
  while (length(x) < n) {
    theta <- runif(n, 0, 2*pi)
    jacobian <- jd_trefoil()
    # Applies the Jacobian scalar value to each value of theta
    jacobian_theta <- sapply(theta, jacobian)
    # Density threshold is the greatest jacobian value in the trefoil knot, and
    # the area of least warping
    density_threshold <- runif(n, 0, jacobian(0))
    # Takes theta values that exceed the density threshold, and throws out the
    # rest
    x <- c(x, theta[jacobian_theta > density_threshold])
  }
  x[1:n]
}

# Jacobian determinant of trefoil
jd_trefoil <- function(){
  function(theta) sqrt(
    17 + 8 * (cos(theta) * cos(2 * theta) - sin(theta) * sin(2 * theta)) +
      9 * (cos(3 * theta))^2
  )
}

Try the tdaunif package in your browser

Any scripts or data that you put into this service are public.

tdaunif documentation built on Sept. 10, 2023, 5:07 p.m.