R/RcppExports.R

Defines functions em_li llike_li log_sum_exp_2 log_sum_exp

Documented in em_li llike_li

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

#' 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 <- function(x) {
    .Call(`_segtest_log_sum_exp`, x)
}

#' 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 <- function(x, y) {
    .Call(`_segtest_log_sum_exp_2`, x, y)
}

#' Objective function for \code{\link{em_li}()}
#'
#' @param B The log-likelihood matrix. Rows are individuals columns are
#'     genotypes.
#' @param lpivec The log prior vector.
#'
#' @author David Gerard
#'
#' @return The log-likelihood of a vector of genotype frequencies when
#'     using genotype likelihoods. This is from Li (2011).
#'
#' @references
#' \itemize{
#'   \item{Li, H. (2011). A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. \emph{Bioinformatics}, 27(21), 2987-2993. \doi{10.1093/bioinformatics/btr509}}
#' }
#'
#' @examples
#' # Simulate some data
#' set.seed(1)
#' gl <- simgl(nvec = c(3, 2, 4, 1, 2))
#' # Log-likelihood at given log-priors
#' prob <- c(0.1, 0.2, 0.4, 0.2, 0.1)
#' lprob <- log(prob)
#' llike_li(B = gl, lpivec = lprob)
#'
#' @export
llike_li <- function(B, lpivec) {
    .Call(`_segtest_llike_li`, B, lpivec)
}

#' EM algorithm from Li (2011)
#'
#' EM algorithm to estimate prior genotype probabilities from genotype
#' likelihoods.
#'
#' @param B Matrix of genotype log-likelihoods. The rows index the individuals
#'     and the columns index the genotypes.
#' @param itermax The maximum number of iterations.
#' @param eps The stopping criteria.
#'
#' @return A vector of log prior probabilities for each genotype.
#'
#' @author David Gerard
#'
#' @references
#' \itemize{
#'   \item{Li, H. (2011). A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. \emph{Bioinformatics}, 27(21), 2987-2993. \doi{10.1093/bioinformatics/btr509}}
#' }
#'
#' @examples
#' # Simulate some data
#' set.seed(1)
#' gl <- simgl(nvec = c(3, 2, 4, 1, 2))
#' # Run em
#' lprob <- em_li(B = gl)
#' # Exponentiate to get probabilities
#' prob <- exp(c(lprob))
#' prob
#'
#' @export
em_li <- function(B, itermax = 100L, eps = 1e-5) {
    .Call(`_segtest_em_li`, B, itermax, eps)
}

Try the segtest package in your browser

Any scripts or data that you put into this service are public.

segtest documentation built on July 1, 2025, 1:07 a.m.