R/dpkb.R

Defines functions dpkb

Documented in dpkb

#' Density function for the Poisson Kernel Based distribution
#' 
#' A function for evaluating the density of the PKB distribution.
#' 
#' @param x vector or matrix of quantiles with length >= 2. If 'x' is a matrix, each row is taken to be a quantile. If 'x' is a vector, it is cast as a matrix with one row.
#' @param mu the location parameter with same length as the quantiles. Normalized to length one.
#' @param rho the concentration parameter, with 0 <= rho < 1.
#' @param logdens logical; if 'TRUE', densities d are given as log(d).
#' @return A vector with the evaluated density values.
#' @export
dpkb <- function(x, mu, rho, logdens = FALSE) {
  # validate input
  if (length(mu) < 2) {
    stop('mu must have length >= 2')
  }
  if (is.vector(x)) {
    x <- matrix(x, ncol = length(x))
  }
  p <- ncol(x)
  if (p < 2) {
    stop('quantiles must have length >= 2')
  }
  if (length(mu) != p) {
    stop('quantiles and mu must have the same length')
  }
  if (rho >= 1 | rho < 0) {
    stop('Input argument rho must be within [0,1)')
  }
  # norm input
  if (sum(abs(mu)) == 0) {
    stop('Input argument mu cannot be a vector of zeros')
  }
  mu <- mu / sqrt(sum(mu^2))
  x <- x / sqrt(rowSums(x^2))

  # calculate log(density)
  logretval <- (log(1 - rho^2) - log(2) - p/2*log(pi) + lgamma(p/2) -
    p/2*log(1 + rho^2 - 2*rho * (x %*% mu)))
  if (logdens) {
    return(logretval)
  } else {
    return(exp(logretval))
  }
}
ahfoss/pkbc documentation built on June 4, 2017, 3:21 a.m.