R/RcppExports.R

In hwep: Hardy-Weinberg Equilibrium in Polyploids

# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

#' Sample gamete counts from full conditional.
#'
#' @param x A vector of genotype counts. x[i] is the number of individuals
#'     with dosage i.
#' @param p A vector of gamete probabilities. p[i] is the proportion of
#'     gametes with dosage i.
#'
#' @author David Gerard
#'
#' @noRd
samp_gametes <- function(x, p) {
    .Call(`_hwep_samp_gametes`, x, p)
}

#' Gibbs sampler under random mating with known genotypes.
#'
#' @param x The vector of genotype counts. x(i) is the number of
#'     individuals that have genotype i.
#' @param alpha Vector of hyperparameters for the gamete frequencies.
#'     Should be length (x.length() - 1) / 2 + 1.
#' @param B The number of sampling iterations.
#' @param T The number of burn-in iterations.
#' @param more A logical. Should we also return posterior draws (\code{TRUE})
#'     or not (\code{FALSE}).
#' @param lg Should we return the log marginal likelihood (true) or not
#'     (false).
#'
#' @return A list with some or all of the following elements
#' \itemize{
#'   \item{\code{mx}: The estimate of the marginal likelihood}
#'   \item{\code{p_tilde}: The value of p used to evaluate the posterior density}.
#'   \item{\code{p}: The samples of the gamete frequencies}
#'   \item{\code{post}: The likelihood times prior evaluated at current samples.}
#'   \item{\code{ptilde_post}: The samples of the full conditionals of p_tilde.}
#' }
#'
#' @author David Gerard
#'
#' @export
gibbs_known <- function(x, alpha, B = 10000L, T = 1000L, more = FALSE, lg = FALSE) {
    .Call(`_hwep_gibbs_known`, x, alpha, B, T, more, lg)
}

#' Calculate marginal likelihood under alternative using genotype likelihoods.
#'
#' Calculates
#' \deqn{
#' \log \prod_i \sum_k l_{ik}q_{ik}
#' }
#' where q = beta / sum(beta). It can return the exponentiated version of this.
#'
#' @param gl genotype log-liklihoods. Rows index individuals, columns
#'     index genotypes.
#' @param beta The concentration parameters.
#' @param lg Should we return the log (true) or the not (false)?
#'
#' @author David Gerard
#'
#' @noRd
#'
plq <- function(gl, beta, lg = FALSE) {
    .Call(`_hwep_plq`, gl, beta, lg)
}

#' Quickly sample from a vector of probabilities.
#'
#' Samples an integer from 0 to probs.length()-1.
#'
#' @param probs The vector of probabilities. Should sum to 1 but it
#'     doesn't check if it does or not, so is very unsafe.
#'
#' @author David Gerard
#'
#' @noRd
sample_int <- function(probs) {
    .Call(`_hwep_sample_int`, probs)
}

#' Sample genotypes from posteriors using genotype likelihoods and genotype priors
#'
#' @param gl The matrix of genotype log-likelihoods. Rows index individuals
#'     and columns index genotypes.
#' @param q The vector of genotype priors (not log-priors).
#'
#' @author David Gerard
#'
#' @noRd
sample_z <- function(gl, q) {
    .Call(`_hwep_sample_z`, gl, q)
}

#' Modify posterior matrix using genotype likelihoods and prior vector.
#'
#' @param postmat The posterior matrix to fill.
#' @param gl The genotype log-likelihoods.
#' @param q The prior genotype probabilities (not logged).
#'
#' @author David Gerard
#'
#' @noRd
mod_postmat <- function(postmat, gl, q) {
    invisible(.Call(`_hwep_mod_postmat`, postmat, gl, q))
}

#' Gibbs sampler under random mating using genotype log-likelihoods.
#'
#' @param gl The matrix of genotype log-likelihoods. The columns index the
#'     dosages and the rows index the individuals. \code{gl[i,j]} is the
#'     genotype log-likelihood for individual i at dosage j. It is assumed
#'     that natural log is used.
#' @param alpha Vector of hyperparameters for the gamete frequencies.
#'     Should be length (x.length() - 1) / 2 + 1.
#' @param B The number of sampling iterations.
#' @param T The number of burn-in iterations.
#' @param more A logical. Should we also return posterior draws (\code{TRUE})
#'     or not (\code{FALSE}).
#' @param lg Should we return the log marginal likelihood (true) or not
#'     (false).
#' @param verbose A logical. Should we print the progress?
#'
#' @return A list with some or all of the following elements
#' \itemize{
#'   \item{\code{mx}: The estimate of the marginal likelihood}
#'   \item{\code{p_tilde}: The value of p used to evaluate the posterior density}.
#'   \item{\code{p}: The samples of the gamete frequencies}
#'   \item{\code{z}: The samples of the individual genotypes}
#'   \item{\code{post}: The samples of the full conditionals of p_tilde.}
#' }
#'
#' @author David Gerard
#'
#' @examples
#' set.seed(1)
#' ploidy <- 8
#'
#' ## Simulate under the null
#' p <- stats::runif(ploidy / 2 + 1)
#' p <- p / sum(p)
#' q <- stats::convolve(p, rev(p), type = "open")
#' nvec <- c(stats::rmultinom(n = 1, size = 100, prob = q))
#' gl <- simgl(nvec)
#'
#' gibbs_gl(gl = gl, alpha = rep(1, ploidy / 2 + 1), lg = TRUE)
#'
#' @export
gibbs_gl <- function(gl, alpha, B = 10000L, T = 1000L, more = FALSE, lg = FALSE, verbose = TRUE) {
    .Call(`_hwep_gibbs_gl`, gl, alpha, B, T, more, lg, verbose)
}

#' Gibbs sampler under the alternative of non-random mating using genotype
#' log-likelihoods.
#'
#' @inheritParams gibbs_gl
#' @param beta The concentration hyperparameter for the genotype frequencies.
#'
#' @return A list with some or all of the following elements
#' \itemize{
#'   \item{\code{mx}: The estimate of the marginal likelihood}
#' }
#'
#' @author David Gerard
#'
#' @examples
#' set.seed(1)
#' ploidy <- 8
#'
#' ## Simulate under the alternative
#' q <- stats::runif(ploidy + 1)
#' q <- q / sum(q)
#' nvec <- c(stats::rmultinom(n = 1, size = 100, prob = q))
#' gl <- simgl(nvec)
#'
#' gibbs_gl_alt(gl = gl, beta = rep(1, ploidy + 1), lg = TRUE)
#'
#' @export
gibbs_gl_alt <- function(gl, beta, B = 10000L, T = 1000L, more = FALSE, lg = FALSE, verbose = TRUE) {
    .Call(`_hwep_gibbs_gl_alt`, gl, beta, B, T, more, lg, verbose)
}

#' Random sample from Dirichlet distribution with n = 1.
#'
#' @param alpha The
#'
#' @author David Gerard
#'
#' @noRd
rdirichlet1 <- function(alpha) {
    .Call(`_hwep_rdirichlet1`, alpha)
}

#' Dirichlet probability density function.
#'
#' @param x The observed vector of proportions. Should sum to 1.
#' @param alpha The concentation parameters, should all be greater than 0.
#' @param lg A logical. Should we return the log pdf or not?
#'
#' @author David Gerard
#'
#' @noRd
ddirichlet <- function(x, alpha, lg = FALSE) {
    .Call(`_hwep_ddirichlet`, x, alpha, lg)
}

#' Multinomial probability mass function
#'
#' @param x The counts.
#' @param p The probabilities.
#' @param lg A logical. Should we return log (true) or not (false).
#'
#' @noRd
#' @author David Gerard
dmultinom_cpp <- function(x, p, lg = FALSE) {
    .Call(`_hwep_dmultinom_cpp`, x, p, lg)
}

#' Log-sum-exponential trick using just two doubles.
#'
#' @param x A double.
#' @param y Another double.
#'
#' @return The log of the sum of the exponential of x and y.
#'
#' @author David Gerard
#'
#' @noRd
log_sum_exp_2_cpp <- function(x, y) {
    .Call(`_hwep_log_sum_exp_2_cpp`, x, y)
}

#' Log-sum-exponential trick.
#'
#' @param x A vector to log-sum-exp.
#'
#' @return The log of the sum of the exponential
#'     of the elements in \code{x}.
#'
#' @author David Gerard
#'
#' @noRd
log_sum_exp_cpp <- function(x) {
    .Call(`_hwep_log_sum_exp_cpp`, x)
}

#' Basic convolution using naive algorithm
#'
#' Same as stats::conv(x, rev(y), type = "open"), but using the naive
#' algorithm. This should be faster for very small vectors (but will
#' be horrible for larger vectors).
#'
#' @param x The first vector
#' @param y The second vector
#'
#' @author David Gerard
#'
#' @noRd
conv_cpp <- function(x, y) {
    .Call(`_hwep_conv_cpp`, x, y)
}