Nothing
#' @title Weight functions
#'
#' @description Marginal density, marginal distribution function,
#' marginal quantile function and random generation for weight functions.
#'
#' @param x,q vector or matrix of quantiles.
#' @param p vector of probabilities.
#' @param n number of observations.
#' @param alpha vector or matrix with concentration parameters for
#' the Dirichlet distribution for a monotonic one.sided or a two.sided
#' weight function.
#' @param alpha1 vector or matrix with concentration parameters for
#' the Dirichlet distribution for the expected direction of non-monotonic
#' one.sided of weight function.
#' @param alpha2 vector or matrix with concentration parameters for
#' the Dirichlet distribution for the unexpected direction of non-monotonic
#' one.sided of weight function.
#' @param omega vector or matrix of fixed probabilities for a one.sided or
#' a two.sided weight function.
#' @param log,log.p logical; if \code{TRUE}, probabilities
#' \code{p} are given as \code{log(p)}.
#' @param lower.tail logical; if \code{TRUE} (default), probabilities
#' are \eqn{P[X \le x]}, otherwise, \eqn{P[X \ge x]}.
#'
#'
#' @examples
#' # draw samples from a two-sided weight function
#' rtwo.sided(10, alpha = c(1, 1))
#'
#' # draw samples from a monotone one-sided weight function
#' rone.sided(10, alpha = c(1, 1, 1))
#'
#' # draw samples from a non-monotone one-sided weight function
#' rone.sided(10, alpha1 = c(1, 1), alpha2 = c(1, 1))
#'
#' @return \code{mdone.sided}, \code{mdtwo.sided}, \code{mdone.sided_fixed},
#' and \code{mdtwo.sided_fixed} give the marginal density,
#' \code{mpone.sided}, \code{mptwo.sided}, \code{mpone.sided_fixed},
#' and \code{mptwo.sided_fixed} give the marginal distribution function,
#' \code{mqone.sided}, \code{mqtwo.sided}, \code{mqone.sided_fixed},
#' and \code{mqtwo.sided_fixed} give the marginal quantile function,
#' and \code{rone.sided}, \code{rtwo.sided}, \code{rone.sided_fixed},
#' and \code{rtwo.sided_fixed} generate random deviates.
#'
#' @export mdone.sided
#' @export mdtwo.sided
#' @export mdone.sided_fixed
#' @export mdtwo.sided_fixed
#' @export rone.sided
#' @export rtwo.sided
#' @export rone.sided_fixed
#' @export rtwo.sided_fixed
#' @export mpone.sided
#' @export mptwo.sided
#' @export mpone.sided_fixed
#' @export mptwo.sided_fixed
#' @export mqone.sided
#' @export mqtwo.sided
#' @export mqone.sided_fixed
#' @export mqtwo.sided_fixed
#' @name weightfunctions
NULL
#### wrappers ####
#' @rdname weightfunctions
mdone.sided <- function(x, alpha = NULL, alpha1 = NULL, alpha2 = NULL, log = FALSE){
# common input check
.check_log(log)
.check_x(x, lower = 0, upper = 1)
if(!is.null(alpha)){
lik <- .mdone.sided_monotonic(x, alpha, log)
}else if(!is.null(alpha1) & !is.null(alpha2)){
lik <- .mdone.sided_general(x, alpha1, alpha2, log)
}
return(lik)
}
#' @rdname weightfunctions
mdtwo.sided <- function(x, alpha, log = FALSE){
# common input check
.check_log(log)
.check_x(x, lower = 0, upper = 1)
lik <- mdone.sided(x, alpha = alpha, log = log)
return(lik)
}
#' @rdname weightfunctions
mdone.sided_fixed <- function(x, omega, log = FALSE){
# common input check
.check_log(log)
.check_x(x, lower = 0, upper = 1)
lik <- .mdone.sided_fixed(x, omega, log)
return(lik)
}
#' @rdname weightfunctions
mdtwo.sided_fixed <- function(x, omega, log = FALSE){
# common input check
.check_log(log)
.check_x(x, lower = 0, upper = 1)
lik <- mdone.sided_fixed(x, omega = omega, log = log)
return(lik)
}
#' @rdname weightfunctions
rone.sided <- function(n, alpha = NULL, alpha1 = NULL, alpha2 = NULL){
# common input check
.check_n(n)
if(!is.null(alpha)){
x <- .rone.sided_monotonic(n, alpha)
}else if(!is.null(alpha1) & !is.null(alpha2)){
x <- .rone.sided_general(n, alpha1, alpha2)
}
return(x)
}
#' @rdname weightfunctions
rtwo.sided <- function(n, alpha){
# common input check
.check_n(n)
x <- rone.sided(n, alpha = alpha)
return(x)
}
#' @rdname weightfunctions
rone.sided_fixed <- function(n, omega){
# common input check
.check_n(n)
x <- .rone.sided_fixed(n, omega)
return(x)
}
#' @rdname weightfunctions
rtwo.sided_fixed <- function(n, omega){
# common input check
.check_n(n)
x <- rone.sided_fixed(n, omega = omega)
return(x)
}
#' @rdname weightfunctions
mpone.sided <- function(q, alpha = NULL, alpha1 = NULL, alpha2 = NULL, lower.tail = TRUE, log.p = FALSE){
# common input check
.check_log.p(log.p)
.check_lower.tail(lower.tail)
.check_q(q, lower = 0, upper = 1)
if(!is.null(alpha)){
p <- .mpone.sided_monotonic(q, alpha, lower.tail, log.p)
}else if(!is.null(alpha1) & !is.null(alpha2)){
p <- .mpone.sided_general(q, alpha1, alpha2, lower.tail, log.p)
}
return(p)
}
#' @rdname weightfunctions
mptwo.sided <- function(q, alpha, lower.tail = TRUE, log.p = FALSE){
# common input check
.check_log.p(log.p)
.check_lower.tail(lower.tail)
.check_q(q, lower = 0, upper = 1)
p <- mpone.sided(q, alpha = alpha, lower.tail = lower.tail, log.p = log.p)
return(p)
}
#' @rdname weightfunctions
mpone.sided_fixed <- function(q, omega, lower.tail = TRUE, log.p = FALSE){
# common input check
.check_log.p(log.p)
.check_lower.tail(lower.tail)
.check_q(q, lower = 0, upper = 1)
p <- .mpone.sided_fixed(q, omega, lower.tail, log.p)
return(p)
}
#' @rdname weightfunctions
mptwo.sided_fixed <- function(q, omega, lower.tail = TRUE, log.p = FALSE){
# common input check
.check_log.p(log.p)
.check_lower.tail(lower.tail)
.check_q(q, lower = 0, upper = 1)
p <- mpone.sided_fixed(q, omega = omega, lower.tail = lower.tail, log.p = log.p)
return(p)
}
#' @rdname weightfunctions
mqone.sided <- function(p, alpha = NULL, alpha1 = NULL, alpha2 = NULL, lower.tail = TRUE, log.p = FALSE){
# common input check
.check_log.p(log.p)
.check_lower.tail(lower.tail)
.check_p(p, log.p)
if(!is.null(alpha)){
q <- .mqone.sided_monotonic(p, alpha, lower.tail, log.p)
}else if(!is.null(alpha1) & !is.null(alpha2)){
q <- .mqone.sided_general(p, alpha1, alpha2, lower.tail, log.p)
}
return(q)
}
#' @rdname weightfunctions
mqtwo.sided <- function(p, alpha, lower.tail = TRUE, log.p = FALSE){
# common input check
.check_log.p(log.p)
.check_lower.tail(lower.tail)
.check_p(p, log.p)
q <- mqone.sided(p, alpha = alpha, lower.tail = lower.tail, log.p = log.p)
return(q)
}
#' @rdname weightfunctions
mqone.sided_fixed <- function(p, omega, lower.tail = TRUE, log.p = FALSE){
# common input check
.check_log.p(log.p)
.check_lower.tail(lower.tail)
.check_p(p, log.p)
q <- .mqone.sided_fixed(p, omega, lower.tail, log.p)
return(q)
}
#' @rdname weightfunctions
mqtwo.sided_fixed <- function(p, omega, lower.tail = TRUE, log.p = FALSE){
# common input check
.check_log.p(log.p)
.check_lower.tail(lower.tail)
.check_p(p, log.p)
q <- mqone.sided_fixed(p, omega = omega, lower.tail = lower.tail, log.p = log.p)
return(q)
}
###### helper functions
#### density functions ####
.mdone.sided_general <- function(x, alpha1, alpha2, log){
stop("Not implemented")
# would require product of two beta-distributed variables, since the weights in the opposite direction
# start at the first cutoff modeled by the first alpha parameter
# https://math.stackexchange.com/questions/1073364/product-of-two-beta-distributed-random-variables
# https://www.dm.fct.unl.pt/sites/www.dm.fct.unl.pt/files/preprints/2012/7_12.pdf
# probably solvable with Meijer g-function
# # input check
# .weightfunctions_check_alpha(alpha1, "alpha1")
# .weightfunctions_check_alpha(alpha2, "alpha2")
#
# # transform to matrices for easier manipulation and checks
# if(length(x) == 1){
# x <- rep(x, nrow(alpha1))
# }
# if(!is.matrix(alpha1)){
# alpha1 <- matrix(alpha1, nrow = 1)
# }
# if(!is.matrix(alpha2)){
# alpha2 <- matrix(alpha2, nrow = 1)
# }
#
# if(nrow(alpha1) != nrow(alpha2))
# stop("Non matching dimensions of 'alpha1' and 'alpha2'.")
# if(nrow(alpha1) != length(x) & nrow(alpha1) != 1)
# stop("Non matching dimensions of 'alpha' and 'x'.")
#
# if(nrow(alpha1) != length(x) & nrow(alpha1) == 1){
# alpha1 <- do.call(rbind, lapply(1:length(x), function(i)alpha1))
# alpha2 <- do.call(rbind, lapply(1:length(x), function(i)alpha2))
# }
#
#
# lik1 <- .mdone.sided_monotonic(x, alpha = alpha1, log = log)
# the side in the unexpected direction starts with the first step of alpha1, e.g., something like this:
# lik2 <- .dmone.sided_monotonic(x, alpha = alpha2, log = log) * beta(alpha[,1], t(apply(alpha[,-1], 1, sum)))
#
# return(cbind(lik2, lik1))
}
.mdone.sided_monotonic <- function(x, alpha, log){
# input check
.weightfunctions_check_alpha(alpha, "alpha")
# transform to matrices for easier manipulation and checks
if(!is.matrix(alpha)){
alpha <- matrix(alpha, nrow = 1)
}
if(length(x) == 1){
x <- rep(x, nrow(alpha))
}
if(nrow(alpha) != length(x) & nrow(alpha) != 1)
stop("Non matching dimensions of 'alpha' and 'x'.")
if(nrow(alpha) != length(x) & nrow(alpha) == 1){
alpha <- do.call(rbind, lapply(1:length(x), function(i)alpha))
}
# marginals of sums of Dirichlet distributed variables are beta distributed
alpha_alpha <- t(apply(alpha, 1, cumsum))
alpha_beta <- t(apply(alpha[,ncol(alpha):1, drop = FALSE], 1, cumsum))[,(ncol(alpha)-1):1, drop = FALSE]
lik <- do.call(cbind, lapply(1:(ncol(alpha) - 1), function(i)stats::dbeta(x, shape1 = alpha_alpha[,i], shape2 = alpha_beta[,i], log = log)))
lik <- cbind(lik, dpoint(x, location = 1))
return(lik)
}
.mdone.sided_fixed <- function(x, omega, log){
.weightfunctions_check_omega(omega, "omega")
# transform to matrices for easier manipulation and checks
if(!is.matrix(omega)){
omega <- matrix(omega, nrow = 1)
}
if(length(x) == 1){
x <- rep(x, nrow(omega))
}
if(nrow(omega) != length(x) & nrow(omega) != 1)
stop("Non matching dimensions of 'omega' and 'x'.")
if(nrow(omega) != length(x) & nrow(omega) == 1){
omega <- do.call(rbind, lapply(1:length(x), function(i)omega))
}
lik <- do.call(cbind, lapply(1:ncol(omega), function(i)dpoint(x, location = omega[,i], log = log)))
}
#### random number generators ####
.rone.sided_general <- function(n, alpha1, alpha2){
# input check
.weightfunctions_check_alpha(alpha1, "alpha1")
.weightfunctions_check_alpha(alpha2, "alpha2")
# transform to matrices for easier manipulation and checks
if(!is.matrix(alpha1)){
alpha1 <- matrix(alpha1, nrow = 1)
}
if(!is.matrix(alpha2)){
alpha2 <- matrix(alpha2, nrow = 1)
}
if(nrow(alpha1) != nrow(alpha2))
stop("Non matching dimensions of 'alpha1' and 'alpha2'.")
if(nrow(alpha1) != n & nrow(alpha1) != 1)
stop("Incompatible dimensions of requested number of samples and 'alpha'.")
if(nrow(alpha1) != n & nrow(alpha1) == 1){
alpha1 <- do.call(rbind, lapply(1:n, function(i)alpha1))
alpha2 <- do.call(rbind, lapply(1:n, function(i)alpha2))
}
x1 <- extraDistr::rdirichlet(n, alpha = alpha1)
x2 <- extraDistr::rdirichlet(n, alpha = alpha2) * (1 - x1[,1])
x <- matrix(ncol = ncol(alpha1) + ncol(alpha2) - 1, nrow = n)
x[,ncol(alpha2):ncol(x)] <- t(apply(x1, 1, cumsum))
for(i in 2:ncol(alpha2)){
x[,i-1] = apply(x2[,i:ncol(alpha2), drop = FALSE], 1, sum) + x1[,1]
}
return(x)
}
.rone.sided_monotonic <- function(n, alpha){
# input check
.weightfunctions_check_alpha(alpha, "alpha")
# transform to matrices for easier manipulation and checks
if(!is.matrix(alpha)){
alpha <- matrix(alpha, nrow = 1)
}
if(nrow(alpha) != n & nrow(alpha) != 1)
stop("Incompatible dimensions of requested number of samples and 'alpha'.")
if(nrow(alpha) != n & nrow(alpha) == 1){
alpha <- do.call(rbind, lapply(1:n, function(i)alpha))
}
x <- extraDistr::rdirichlet(n, alpha = alpha)
x <- t(apply(x, 1, cumsum))
return(x)
}
.rone.sided_fixed <- function(n, omega){
.weightfunctions_check_omega(omega, "omega")
# transform to matrices for easier manipulation and checks
if(!is.matrix(omega)){
omega <- matrix(omega, nrow = 1)
}
if(nrow(omega) != n & nrow(omega) != 1)
stop("Incompatible dimensions of requested number of samples and 'omega'.")
if(nrow(omega) != n & nrow(omega) == 1){
omega <- do.call(rbind, lapply(1:n, function(i)omega))
}
x <- omega
return(x)
}
#### marginal distribution functions ####
.mpone.sided_general <- function(q, alpha1, alpha2, lower.tail, log.p){
stop("Not implemented")
# # input check
# .weightfunctions_check_alpha(alpha1, "alpha1")
# .weightfunctions_check_alpha(alpha2, "alpha2")
#
#
# # transform to matrices for easier manipulation and checks
# if(!is.matrix(alpha1)){
# alpha1 <- matrix(alpha1, nrow = 1)
# }
# if(!is.matrix(alpha2)){
# alpha2 <- matrix(alpha2, nrow = 1)
# }
# if(length(q) == 1){
# q <- rep(q, nrow(alpha1))
# }
#
# if(nrow(alpha1) != nrow(alpha2))
# stop("Non matching dimensions of 'alpha1' and 'alpha2'.")
# if(nrow(alpha1) != length(q) & nrow(alpha1) != 1)
# stop("Non matching dimensions of 'alpha' and 'q'.")
}
.mpone.sided_monotonic <- function(q, alpha, lower.tail, log.p){
# input check
.weightfunctions_check_alpha(alpha, "alpha")
# transform to matrices for easier manipulation and checks
if(!is.matrix(alpha)){
alpha <- matrix(alpha, nrow = 1)
}
if(length(q) == 1){
q <- rep(q, nrow(alpha))
}
if(nrow(alpha) != length(q) & nrow(alpha) != 1)
stop("Non matching dimensions of 'alpha' and 'q'.")
if(nrow(alpha) != length(q) & nrow(alpha) == 1){
alpha <- do.call(rbind, lapply(1:length(q), function(i)alpha))
}
# marginals of sums of Dirichlet distributed variables are beta distributed
alpha_alpha <- t(apply(alpha, 1, cumsum))
alpha_beta <- t(apply(alpha[,ncol(alpha):1, drop = FALSE], 1, cumsum))[,(ncol(alpha)-1):1, drop = FALSE]
p <- do.call(cbind, lapply(1:(ncol(alpha) - 1), function(i)stats::pbeta(q, shape1 = alpha_alpha[,i], shape2 = alpha_beta[,i], log.p = log.p, lower.tail = lower.tail)))
p <- cbind(p, ppoint(q, location = 1, lower.tail = lower.tail, log.p = log.p))
return(p)
}
.mpone.sided_fixed <- function(q, omega, lower.tail, log.p){
.weightfunctions_check_omega(omega, "omega")
# transform to matrices for easier manipulation and checks
if(!is.matrix(omega)){
omega <- matrix(omega, nrow = 1)
}
if(length(q) == 1){
q <- rep(q, nrow(omega))
}
if(nrow(omega) != length(q) & nrow(omega) != 1)
stop("Non matching dimensions of 'omega' and 'q'.")
if(nrow(omega) != length(q) & nrow(omega) == 1){
omega <- do.call(rbind, lapply(1:length(q), function(i)omega))
}
p <- do.call(cbind, lapply(1:ncol(omega), function(i)ppoint(q, location = omega[,i], lower.tail = lower.tail, log.p = log.p)))
return(p)
}
#### marginal quantile functions ####
.mqone.sided_general <- function(p, alpha1, alpha2, lower.tail, log.p){
stop("Not implemented")
# # input check
# .weightfunctions_check_alpha(alpha1, "alpha1")
# .weightfunctions_check_alpha(alpha2, "alpha2")
#
#
# # transform to matrices for easier manipulation and checks
# if(!is.matrix(alpha1)){
# alpha1 <- matrix(alpha1, nrow = 1)
# }
# if(!is.matrix(alpha2)){
# alpha2 <- matrix(alpha2, nrow = 1)
# }
# if(length(p) == 1){
# p <- rep(p, nrow(alpha1))
# }
#
# if(nrow(alpha1) != nrow(alpha2))
# stop("Non matching dimensions of 'alpha1' and 'alpha2'.")
# if(nrow(alpha1) != length(p) & nrow(alpha1) != 1)
# stop("Non matching dimensions of 'alpha' and 'p'.")
#
# if(nrow(alpha1) != length(p) & nrow(alpha1) == 1){
# alpha1 <- do.call(rbind, lapply(1:length(p), function(i)alpha1))
# alpha2 <- do.call(rbind, lapply(1:length(p), function(i)alpha2))
# }
}
.mqone.sided_monotonic <- function(p, alpha, lower.tail, log.p){
# input check
.weightfunctions_check_alpha(alpha, "alpha")
# transform to matrices for easier manipulation and checks
if(!is.matrix(alpha)){
alpha <- matrix(alpha, nrow = 1)
}
if(length(q) == 1){
p <- rep(p, nrow(alpha))
}
if(nrow(alpha) != length(p) & nrow(alpha) != 1)
stop("Non matching dimensions of 'alpha' and 'p'.")
if(nrow(alpha) != length(p) & nrow(alpha) == 1){
alpha <- do.call(rbind, lapply(1:length(p), function(i)alpha))
}
# marginals of sums of Dirichlet distributed variables are beta distributed
alpha_alpha <- t(apply(alpha, 1, cumsum))
alpha_beta <- t(apply(alpha[,ncol(alpha):1, drop = FALSE], 1, cumsum))[,(ncol(alpha)-1):1, drop = FALSE]
q <- do.call(cbind, lapply(1:(ncol(alpha) - 1), function(i)stats::qbeta(p, shape1 = alpha_alpha[,i], shape2 = alpha_beta[,i], lower.tail = lower.tail, log.p = log.p)))
q <- cbind(q, qpoint(p, location = 1, lower.tail = lower.tail, log.p = log.p))
# make sure that the range is correct qpoint operates on unrestricted range
q <- ifelse(q == -Inf, 0, q)
q <- ifelse(q == Inf, 1, q)
return(q)
}
.mqone.sided_fixed <- function(p, omega, lower.tail, log.p){
.weightfunctions_check_omega(omega, "omega")
# transform to matrices for easier manipulation and checks
if(!is.matrix(omega)){
omega <- matrix(omega, nrow = 1)
}
if(length(q) == 1){
p <- rep(p, nrow(omega))
}
if(nrow(omega) != length(p) & nrow(omega) != 1)
stop("Non matching dimensions of 'omega' and 'p'.")
if(nrow(omega) != length(p) & nrow(omega) == 1){
omega <- do.call(rbind, lapply(1:length(p), function(i)omega))
}
q <- do.call(cbind, lapply(1:ncol(omega), function(i)qpoint(p, location = omega[,i], lower.tail = lower.tail, log.p = log.p)))
# make sure that the range is correct
q <- ifelse(q == -Inf, 0, q)
q <- ifelse(q == Inf, 1, q)
return(q)
}
### helper functions
.weightfunctions_check_alpha <- function(alpha, name = "alpha"){
if(!is.numeric(alpha) | !(is.vector(alpha) | is.matrix(alpha)))
stop(paste0("'", name, "' must be a numeric vector or a matrix."))
if(is.vector(alpha))if(length(alpha) < 2)
stop(paste0("'", name, "' must be a vector of length at least 2."))
if(is.matrix(alpha))if(ncol(alpha) < 2)
stop(paste0("'", name, "' must be a matrix with at least 2 columns."))
if(!all(alpha > 0))
stop(paste0("'", name, "' must be positive."))
}
.weightfunctions_check_omega <- function(omega, name = "omega"){
if(!is.numeric(omega) | !(is.vector(omega) | is.matrix(omega)))
stop(paste0("'", name, "' must be a numeric vector or a matrix."))
if(is.vector(omega))if(length(omega) < 2)
stop(paste0("'", name, "' must be a vector of length at least 2."))
if(is.matrix(omega))if(ncol(omega) < 2)
stop(paste0("'", name, "' must be a matrix with at least 2 columns."))
if(!all(omega >= 0 & omega <= 1))
stop(paste0("'", name, "' must be between 0 and 1."))
}
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.