Nothing
#' @title The Elliptic Distribution
#'
#' @description
#' Density, distribution function, quantile function,
#' and random generation for the univariate elliptic distribution.
#'
#' @param p numeric vector of probabilities.
#' @param x numeric vector of quantiles.
#' @param q numeric vector of quantiles.
#' @param n number of observations.
#' @param object an object of \code{\link{ecd-class}}.
#' To achieve high performance for \code{qec} and \code{rec},
#' it should be created with \code{with.quantile=TRUE}.
#' @param debug logical, whether to print debug message, default is \code{FALSE}.
#'
#' @return \code{dec} gives the density,
#' \code{pec} gives the distribution function,
#' \code{qec} gives the quantile function,
#' \code{rec} generates random deviates.
#'
#' @keywords distribution
#'
#' @include ecd-class.R
#'
#' @export dec
#' @export pec
#' @export qec
#' @export rec
#'
#' @importFrom stats predict runif
#'
#' @author Stephen H. Lihn
#'
#' @examples
#' d <- ecd(with.quantile=TRUE)
#' x <- seq(-20, 20, by=5)
#' dec(x,d)
#' pec(x,d)
#' p <- c(0.0001, 0.001, 0.01, 0.99, 0.999, 0.9999)
#' qec(p,d)
#' rec(100,d)
### <======================================================================>
"dec" <- function(x, object=ecd()) {
ecd.pdf(object, x)
}
#' @rdname dec
"pec" <- function(q, object = ecd()) {
if (length(object@stats)==0) {
stop("stats is not computed in ecd object")
}
N <- 4
mean <- object@stats$mean
dv <- object@stats$stdev
xt1 <- mean - N*dv # left tail
xt2 <- mean + N*dv # right tail
# (a) switch piece.wise off for tail region
# (b) switch to ccdf for q > mean
f1 <- function(q) ecd.cdf(object, q, piece.wise=TRUE)
f2 <- function(q) { 1-ecd.ccdf(object, q, piece.wise=TRUE) }
f3 <- function(q) ecd.cdf(object, q, piece.wise=FALSE)
f4 <- function(q) { 1-ecd.ccdf(object, q, piece.wise=FALSE) }
Nq <- length(q)
p <- ecd.mpnum(object, rep(NaN, Nq))
insert <- function(i,v) { p[i] <<- v }
i1 <- which(q <= mean & q > xt1)
i2 <- which(q > mean & q < xt2)
i3 <- which(q <= xt1)
i4 <- which(q >= xt2)
if (length(i1)>0) mapply(insert, i1, f1(q[i1]))
if (length(i2)>0) mapply(insert, i2, f2(q[i2]))
if (length(i3)>0) mapply(insert, i3, f3(q[i3]))
if (length(i4)>0) mapply(insert, i4, f4(q[i4]))
ecd.mpnum(object, p)
}
#' @rdname dec
"qec" <- function(p, object = ecd(with.quantile=TRUE), debug=FALSE) {
if (! ecd.has_quantile(object)) {
object <- quantilize(object, show.warning=TRUE)
}
dq <- object@quantile
p <- ecd.mp2f(p) # lm can't handle mpfr, so revert to numeric
in_range <- (p>=0 & p<=1)
if (! all(in_range)) {
stop("Probability input out of range (0,1)!")
}
between_seg <- function(idx, c.from, c.to, q) {
if (q >= c.from & q <= c.to) idx else 0
}
# translate q to segment index
get_idx <- function(q) {
max(mapply(between_seg,
seq(1, dq@N_seg),
dq@cseg.from, dq@cseg.to,
MoreArgs = list(q=q),
SIMPLIFY = TRUE))
}
predict_quantile <- function(q) {
idx <- get_idx(q)
if (debug) print(paste("idx", idx, "of", dq@N_seg, "for q=", q))
if (idx > 0) {
if (debug) {
print(paste("xseg from", dq@xseg.from[idx], " to ", dq@xseg.to[idx]))
print(paste("cseg from", dq@cseg.from[idx], " to ", dq@cseg.to[idx]))
}
fit <- dq@cdf.fit[[idx]]
if (debug) print(fit)
xp <- predict(fit, data.frame(cdf=c(q)))
} else {
# tail region
if (q <= dq@cseg.min & q > 0) {
# left tail
if (debug) {
print(paste("left tail, log cdf for q=", q))
print(dq@fit.left)
}
xp <- predict(dq@fit.left, data.frame(lcdf=c(log(q))))
} else if (q >= dq@cseg.max & q < 1) {
# right tail
if (debug) {
print(paste("right tail, log ccdf for 1-q, q=", q))
print(dq@fit.right)
}
xp <- predict(dq@fit.right, data.frame(lccdf=c(log(1-q))))
}
else if (q==0) { # edge case
if (debug) print(paste("edge case for q=", q))
xp <- -Inf
} else if (q==1) {
if (debug) print(paste("edge case for q=", q))
xp <- Inf
} else {
stop(paste("Failed to locate quantile region. q=",q))
}
}
xp
}
# quantile specified, give me x
if (debug & length(p)==1) {
return(ecd.mpnum(object, predict_quantile(p)))
}
rs <- simplify2array(parallel::mclapply(p, predict_quantile))
ecd.mpnum(object, rs)
}
#' @rdname dec
"rec" <- function(n, object = ecd(with.quantile=TRUE)) {
qec(runif(n), object)
}
### <---------------------------------------------------------------------->
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.