R/detleafs.R

In detpack: Density Estimation and Random Number Generation with Distribution Element Trees

#' Extract Leaf Elements from Distribution Element Tree
#'
#' The function \code{det.leafs} extracts the distribution elements at the branch ends of a DET generated by the function \code{\link{det.construct}}.
#'
#' @param det distribution element tree object resulting from \code{\link{det.construct}}.
#'
#' @return A list of vectors containing the leaf element data is returned: \code{p} probability density, \code{theta} element parameters, \code{lb} lower bound, \code{size} of element, \code{div} divisions or splits along dimensions leading to final element.
#' @export
#'
#' @examples
#' require(stats); require(graphics)
#' # generate DET based on bi-variate Gaussian data
#' n <- 1e4; x <- rnorm(n)
#' x <- matrix(c(x, x+rnorm(n,0,0.2)), nrow = 2, byrow = TRUE)
#' det <- det.construct(x)
#' # plot data and element pattern
#' leafs <- det.leafs(det)
#' plot(t(x), type = "p", pch = ".", asp = 1)
#' for (k in 1:length(leafs$p)) {
#'    p <- leafs$lb[,k] # element corner point
#'    w <- leafs$size[,k] # element size
#'    elem <- rbind(c(p[1],p[1]+w[1],p[1]+w[1],p[1],p[1]),
#'                  c(p[2],p[2],p[2]+w[2],p[2]+w[2],p[2])) # element rectangle
#'    elem <- t(det$A) %*% elem + det$mu %*% t(rep(1,5)) # pre-white transform
#'    lines(elem[1,],elem[2,]) # draw element
#' }
det.leafs <- function(det) {
   d <- length(det$lb) # number of dimensions
   indl <- seq(nrow(det$tree))[is.na(det$tree[,2])] # indices of leafs in tree
   m <- length(indl) # number of leaf elements
   # initialize leaf property fields
   p <- rep(NA,m) # de probability density
   theta <- rep(list(NA),m) # de parameters
   lb <- matrix(rep(NA,d*m), nrow = d) # de lower bound
   size <- lb # de size
   div <- lb # de subdivisions
   # extract leaf data
   for (k in 1:m) { # loop over leaf elements
      dek <- det.de(det, indl[k])
      p[k] <- dek$p; theta[[k]] <- dek$theta
      lb[,k] <- dek$lb; size[,k] <- dek$size; div[,k] <- dek$div
   }
   return(list(p = p, theta = theta, lb = lb, size = size, div = div))
}

#' Extract Distribution Element Characteristics
#'
#' The function \code{det.de} extracts the distribution element with index \code{ind} from a distribution element tree (DET) generated by the function \code{\link{det.construct}}.
#'
#' @param det distribution element tree object resulting from \code{\link{det.construct}}.
#' @param ind index of element to extract from \code{det}.
#'
#' @return A list with the element characteristics is returned: \code{p} probability density, \code{theta} element parameters, \code{lb} lower bound, \code{size} of element, \code{div} divisions or splits along dimensions leading to final element.
#' @export
det.de <- function(det, ind) {
   d <- length(det$lb) # number of dimensions
   # initialize de property fields
   lb <- rep(0,d) # lower bound
   size <- rep(1,d) # size
   div <- rep(0,d) # subdivisions or splits
   # extract de data
   p <- det$p[ind] # probability
   theta <- det$theta[[ind]] # parameters
   # determine de center, size, and subdivisions
   while (!is.na(det$tree[ind,1])) { # move from de to root de
      indp <- det$tree[ind,1] # move one level down to parent
      dimens <- det$sd[indp]; pos <- det$sp[indp] # split dimension & position
      div[dimens] <- div[dimens]+1
      if (det$tree[indp,2] == ind) { # first child
         size[dimens] <- size[dimens]*pos
         lb[dimens] <- lb[dimens]*pos
      } else { # second child
         size[dimens] <- size[dimens]*(1-pos)
         lb[dimens] <- lb[dimens]*(1-pos) + pos
      }
      ind <- indp # move one level down to parent
   }
   # rescale
   size <- size * (det$ub - det$lb)
   lb <- det$lb + lb * (det$ub - det$lb)
   p <- p/prod(det$ub - det$lb)
   return(list(p = p, theta = theta, lb = lb, size = size, div = div))
}