R/RcppExports.R

Defines functions bbpost_double bbpost_tot get_out_prop dbeta_dprop dbeta_dh dh_dr dd_ds dbeta_dr dxi_df df_deps deps_dell dbeta_dl dxi_dd dbeta_ds dbeta_dr_ell grad_offspring_mat grad_offspring grad_offspring_weights dbeta_deps dbeta_dd dh_dtau dbeta_dtau grad_offspring_mat_original grad_offspring_original grad_offspring_weights_original obj_offspring_vec obj_offspring obj_offspring_reparam obj_offspring_weights obj_offspring_weights_reparam outlier_obj dbeta_dtau_withxi outlier_grad obj_parent obj_parent_reparam grad_parent grad_parent_reparam out_grad_parent get_parent_outprop pbias pbias_double dbetabinom_cpp dbetabinom_mu_rho_cpp dbetabinom_mu_rho_cpp_double dhyper_cpp get_prob_geno get_q_array_cpp logsumexp get_pvec colSums_cpp expit

Documented in bbpost_double bbpost_tot colSums_cpp dbetabinom_cpp dbetabinom_mu_rho_cpp dbetabinom_mu_rho_cpp_double dbeta_dd dbeta_deps dbeta_dh dbeta_dl dbeta_dprop dbeta_dr dbeta_dr_ell dbeta_ds dbeta_dtau dbeta_dtau_withxi dd_ds deps_dell df_deps dh_dr dh_dtau dhyper_cpp dxi_dd dxi_df expit get_out_prop get_parent_outprop get_prob_geno get_pvec get_q_array_cpp grad_offspring grad_offspring_mat grad_offspring_mat_original grad_offspring_original grad_offspring_weights grad_offspring_weights_original grad_parent grad_parent_reparam logsumexp obj_offspring obj_offspring_reparam obj_offspring_vec obj_offspring_weights obj_offspring_weights_reparam obj_parent obj_parent_reparam out_grad_parent outlier_grad outlier_obj pbias pbias_double

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

#' Posterior inference under the beta-binomial model with outliers. If \code{od_param = 0}, then
#' this reduces to the binomial model with outliers.
#'
#' @param x The observed number of reference alleles.
#' @param n The total number of reads.
#' @param ploidy The ploidy of the species
#' @param prob_geno The distribution of genotypes.
#' @param seq_error The sequencing error rate.
#' @param bias_val The bias parameter. A value of 1 indicates no bias.
#' @param od_param The over-dispersion parameter. Should be in [0, 1). A value of 0 reduces to the
#'     binomial model.
#' @param outlier A logical. Should we include the outlier model (\code{TRUE}) or not (\code{FALSE})?
#' @param out_prop The proportion of observations that are outliers.
#' @param out_mean The mean of outlier distribution.
#' @param out_disp The overdispersion parameter of the outlier distribution.
#'
#' @author David Gerard
#'
#' @export
#'
bbpost_double <- function(x, n, ploidy, prob_geno, seq_error, bias_val, od_param, outlier = FALSE, out_prop = 0.1, out_mean = 0.5, out_disp = 1/3) {
    .Call('_updogAlpha_bbpost_double', PACKAGE = 'updogAlpha', x, n, ploidy, prob_geno, seq_error, bias_val, od_param, outlier, out_prop, out_mean, out_disp)
}

#' Posterior inference for each individual.
#'
#' @inheritParams bbpost_double
#' @param ocounts The vector of counts of the reference allele of the offspring.
#' @param osize The vector of total reads of the offspring.
#'
#' @author David Gerard
#'
#' @export
#'
bbpost_tot <- function(ocounts, osize, ploidy, prob_geno, seq_error, bias_val, od_param, outlier = FALSE, out_prop = 0.1, out_mean = 0.5, out_disp = 1/3) {
    .Call('_updogAlpha_bbpost_tot', PACKAGE = 'updogAlpha', ocounts, osize, ploidy, prob_geno, seq_error, bias_val, od_param, outlier, out_prop, out_mean, out_disp)
}

#' E step in EM algorithm for a single individual.
#'
#' @inheritParams obj_offspring
#' @inheritParams dbeta_deps
#'
#'
#' @author David Gerard
#'
get_out_prop <- function(ocounts, osize, ploidy, prob_geno, d, eps, tau, out_prop, out_mean, out_disp) {
    .Call('_updogAlpha_get_out_prop', PACKAGE = 'updogAlpha', ocounts, osize, ploidy, prob_geno, d, eps, tau, out_prop, out_mean, out_disp)
}

#' Derivative of mean of beta density.
#'
#' This calculates d/dxi BB(x|n, xi, tau).
#'
#' @param x The observed counts of reference allele.
#' @param n The total number of observed counts.
#' @param xi The mean proportion.
#' @param tau The overdispersion parameter.
#'
#' @author David Gerard
#'
dbeta_dprop <- function(x, n, xi, tau) {
    .Call('_updogAlpha_dbeta_dprop', PACKAGE = 'updogAlpha', x, n, xi, tau)
}

#' Derivative of overdispersion parameter of beta density.
#'
#' This calculates d/dh of BB(x|n, xi, h), where h = (1 - tau) / tau, so tau = 1 / (h + 1).
#'
#' @inheritParams dbeta_dprop
#' @param h A parameterization of the overdispersion parameter.
#'     This is just constrained to be greater than 0.
#'
#' @author David Gerard
#'
dbeta_dh <- function(x, n, xi, h) {
    .Call('_updogAlpha_dbeta_dh', PACKAGE = 'updogAlpha', x, n, xi, h)
}

#' Just a wrapper for \code{std::exp}.
#'
#' @author David Gerard
#'
#' @param r A double.
#'
#'
dh_dr <- function(r) {
    .Call('_updogAlpha_dh_dr', PACKAGE = 'updogAlpha', r)
}

#' Just a wrapper for \code{std::exp}.
#'
#' @author David Gerard
#'
#' @param s A double.
#'
#'
dd_ds <- function(s) {
    .Call('_updogAlpha_dd_ds', PACKAGE = 'updogAlpha', s)
}

#' Derivative of beta(x|n, xi, r), where r = log(h) from \code{\link{dbeta_dh}}.
#'
#' @inheritParams dbeta_dh
#' @param r We have \code{r = log(h)} from double dbeta_dr().
#'
#' @author David Gerard
#'
dbeta_dr <- function(x, n, xi, r) {
    .Call('_updogAlpha_dbeta_dr', PACKAGE = 'updogAlpha', x, n, xi, r)
}

#' Returns derivative to f / (d * (1 - f) + f)
#'
#' @param f The sequencing-error adjusted probability of the refernence allele.
#' @param d The bias term. Must be greater than 0.
#'
#' @author David Gerard
#'
dxi_df <- function(d, f) {
    .Call('_updogAlpha_dxi_df', PACKAGE = 'updogAlpha', d, f)
}

#' Returns derivative of p(1 - eps) + (1 - p) * eps
#'
#' @param p The proportion of genomes whose allele is A.
#' @param eps The sequencing error rate.
#'
#'
#' @author David Gerard
#'
df_deps <- function(eps, p) {
    .Call('_updogAlpha_df_deps', PACKAGE = 'updogAlpha', eps, p)
}

#' Derivative of exp(ell) / (1 + exp(ell))
#'
#' exp(ell) / (1 + exp(ell)) is the sequencing error rate and ell is a parameterization to
#' that is unconstrained.
#'
#' @param ell A double.
#'
#' @author David Gerard
#'
deps_dell <- function(ell) {
    .Call('_updogAlpha_deps_dell', PACKAGE = 'updogAlpha', ell)
}

#' Derivative of beta density w.r.t. unconstrained parameterization of sequencing error.
#'
#' I use the chain rule here. \code{\link{dbeta_dprop}} * \code{\link{dxi_df}} *
#' \code{\link{df_deps}} * \code{\link{deps_dell}}.
#'
#' @param x The number of counts of reference allele.
#' @param n The number of counts of reads.
#' @param d The sequencing bias parameter.
#' @param ell Expit of eps. We have eps = exp(ell) / (1 + exp(ell)), or ell = log(eps / (1 - eps))
#' @param p The proportion of genome that is the reference allele.
#' @param h the overdisperion parameter. tau = 1 / (h + 1) or h = (1 - tau) / tau
#'
#'
#' @author David Gerard
#'
dbeta_dl <- function(x, n, d, ell, p, h) {
    .Call('_updogAlpha_dbeta_dl', PACKAGE = 'updogAlpha', x, n, d, ell, p, h)
}

#' Derivative w.r.t. d of xi(d, f) = f / (d * (1 - f) + f)
#'
#' @inheritParams dxi_df
#'
#' @author David Gerard
#'
#'
dxi_dd <- function(d, f) {
    .Call('_updogAlpha_dxi_dd', PACKAGE = 'updogAlpha', d, f)
}

#' Derivative of betabinomial density w.r.t. bias parameter.
#'
#' Uses chain rule with \code{\link{dbeta_dprop}} * \code{\link{dxi_dd}} * \code{\link{dd_ds}}.
#'
#' @inheritParams dbeta_dl
#' @param s We have \code{s = exp(d)}, where \code{d} is the bias parameter.
#'
#' @author David Gerard
#'
dbeta_ds <- function(x, n, s, ell, p, h) {
    .Call('_updogAlpha_dbeta_ds', PACKAGE = 'updogAlpha', x, n, s, ell, p, h)
}

#' Same as \code{\link{dbeta_dh}}, but with same inputs as \code{\link{dbeta_ds}}
#' and \code{\link{dbeta_dl}}.
#'
#' @inheritParams dbeta_dl
#' @param r We have \code{tau = 1 / (1 + exp(r))}
#'
#' @author David Gerard
#'
dbeta_dr_ell <- function(x, n, d, ell, p, r) {
    .Call('_updogAlpha_dbeta_dr_ell', PACKAGE = 'updogAlpha', x, n, d, ell, p, r)
}

#' Gradient of \code{\link{obj_offspring_reparam}} for each individual.
#'
#' @inheritParams obj_offspring
#' @param s We have \code{s = exp(d)}, where \code{exp(d)} is the bias term.
#' @param ell The logit of the sequencing error rate
#' @param r We have \code{r = log((1 - tau) / tau)}
#'
#'
#' @author David Gerard
#'
grad_offspring_mat <- function(ocounts, osize, ploidy, prob_geno, s, ell, r) {
    .Call('_updogAlpha_grad_offspring_mat', PACKAGE = 'updogAlpha', ocounts, osize, ploidy, prob_geno, s, ell, r)
}

#' Gradient of \code{\link{obj_offspring_reparam}}.
#'
#' @inheritParams obj_offspring
#' @inheritParams grad_offspring_mat
#'
#' @author David Gerard
#'
#' @export
#'
grad_offspring <- function(ocounts, osize, ploidy, prob_geno, s, ell, r) {
    .Call('_updogAlpha_grad_offspring', PACKAGE = 'updogAlpha', ocounts, osize, ploidy, prob_geno, s, ell, r)
}

#' Gradient of \code{\link{obj_offspring_weights_reparam}}
#'
#' @inheritParams grad_offspring
#' @param weight_vec A vector of weights between 0 and 1 (do not need to add up to 1).
#'
grad_offspring_weights <- function(ocounts, osize, weight_vec, ploidy, prob_geno, s, ell, r) {
    .Call('_updogAlpha_grad_offspring_weights', PACKAGE = 'updogAlpha', ocounts, osize, weight_vec, ploidy, prob_geno, s, ell, r)
}

#' Derivative of beta density w.r.t. sequencing error rate.
#'
#' I use the chain rule here. \code{\link{dbeta_dprop}} * \code{\link{dxi_df}} *
#' \code{\link{df_deps}}.
#'
#' @param x The number of counts of reference allele.
#' @param n The number of counts of reads.
#' @param d The sequencing bias parameter.
#' @param eps The sequencing error rate
#' @param p The proportion of genome that is the reference allele.
#' @param tau The overdisperion parameter.
#'
#'
#' @author David Gerard
#'
dbeta_deps <- function(x, n, d, eps, p, tau) {
    .Call('_updogAlpha_dbeta_deps', PACKAGE = 'updogAlpha', x, n, d, eps, p, tau)
}

#' Derivative of betabinomial density w.r.t. original bias parameter.
#'
#' Uses chain rule with \code{\link{dbeta_dprop}} * \code{\link{dxi_dd}}.
#'
#' @inheritParams dbeta_deps
#'
#' @author David Gerard
#'
dbeta_dd <- function(x, n, d, eps, p, tau) {
    .Call('_updogAlpha_dbeta_dd', PACKAGE = 'updogAlpha', x, n, d, eps, p, tau)
}

#' Derivative of h(tau) = (1 - tau) / tau.
#'
#' Just returns -1 / tau^2
#'
#' @param tau A double. The overdispersion parameter
#'
#' @author David Gerard
#'
dh_dtau <- function(tau) {
    .Call('_updogAlpha_dh_dtau', PACKAGE = 'updogAlpha', tau)
}

#' Derivative of beta(x|n, xi, tau).
#'
#' @inheritParams dbeta_deps
#'
#' @author David Gerard
#'
dbeta_dtau <- function(x, n, d, eps, p, tau) {
    .Call('_updogAlpha_dbeta_dtau', PACKAGE = 'updogAlpha', x, n, d, eps, p, tau)
}

#' Gradient of \code{\link{obj_offspring}} for each individual. This is in the
#' original parameterization.
#'
#' @inheritParams obj_offspring
#' @inheritParams dbeta_deps
#'
#'
#' @author David Gerard
#'
grad_offspring_mat_original <- function(ocounts, osize, ploidy, prob_geno, d, eps, tau) {
    .Call('_updogAlpha_grad_offspring_mat_original', PACKAGE = 'updogAlpha', ocounts, osize, ploidy, prob_geno, d, eps, tau)
}

#' Gradient of \code{\link{obj_offspring}} using the original parameterization.
#'
#' @inheritParams obj_offspring
#' @inheritParams dbeta_deps
#'
#' @return A NumericVector of length three with the partial derivatives of
#'     \code{d}, \code{eps}, and \code{tau}, in that order.
#'
#' @author David Gerard
#'
#' @export
#'
grad_offspring_original <- function(ocounts, osize, ploidy, prob_geno, d, eps, tau) {
    .Call('_updogAlpha_grad_offspring_original', PACKAGE = 'updogAlpha', ocounts, osize, ploidy, prob_geno, d, eps, tau)
}

#' Gradient of \code{\link{obj_offspring_weights}} using original parameterization
#'
#' @inheritParams obj_offspring
#' @inheritParams dbeta_deps
#' @param weight_vec A vector of weights between 0 and 1 (do not need to add up to 1).
#'
#' @return A NumericVector of length three with the partial derivatives of
#'     \code{d}, \code{eps}, and \code{tau}, in that order.
#'
grad_offspring_weights_original <- function(ocounts, osize, weight_vec, ploidy, prob_geno, d, eps, tau) {
    .Call('_updogAlpha_grad_offspring_weights_original', PACKAGE = 'updogAlpha', ocounts, osize, weight_vec, ploidy, prob_geno, d, eps, tau)
}

#' Vector of objective functions for offspring.
#'
#' @param ocounts The observed counts of the refernce
#'     allele for each individual.
#' @param osize The observed number of reads for each
#'     individuals.
#' @param ploidy An integer. The ploidy of the species. This is assumed
#'     to be the same for each individual.
#' @param prob_geno The allele frequencies of the genotypes. See \code{\link{get_prob_geno}}.
#' @param bias_val The bias parameter. A value of 1 means there is no bias
#'     toward one allele or the other. A value less than one indicates a bias
#'     toward the reference allele. A value greater than one indicates a bias
#'     toward the non-reference allele.
#' @param seq_error The sequencing error rate.
#' @param od_param The overdispersion parameter in the beta-binomial model
#'     for the OK counts. When this is zero, this resorts to the binomial
#'     model for counts.
#' @param outlier A logical. Should we include an outlier model (\code{TRUE})
#'     or not (\code{FALSE})? Defaults to \code{FALSE}.
#' @param out_prop The proportion of points that are outliers. Defaults
#'     (quite arbitrarily) to \code{0.01}.
#' @param out_mean The mean of beta-binomial for the outlier distribution.
#'     Defaults to \code{0.5}.
#' @param out_disp The overdispersion parameter of the outlier distribution.
#'     Defaults to \code{1/3}, which corresponds to a uniform distribution
#'     for the underlying beta when \code{out_mean = 0.5}.
#'
#' @author David Gerard
#'
#'
#' @seealso \code{\link{up_bb_obj}}.
#'
obj_offspring_vec <- function(ocounts, osize, ploidy, prob_geno, bias_val = 1, seq_error = 0, od_param = 0, outlier = FALSE, out_prop = 0.01, out_mean = 0.5, out_disp = 1.0 / 3.0) {
    .Call('_updogAlpha_obj_offspring_vec', PACKAGE = 'updogAlpha', ocounts, osize, ploidy, prob_geno, bias_val, seq_error, od_param, outlier, out_prop, out_mean, out_disp)
}

#' Objective function for the offspring.
#'
#' @inheritParams obj_offspring_vec
#'
#' @author David Gerard
#'
#' @export
#'
obj_offspring <- function(ocounts, osize, ploidy, prob_geno, bias_val = 1, seq_error = 0, od_param = 0, outlier = FALSE, out_prop = 0.01, out_mean = 0.5, out_disp = 1.0 / 3.0) {
    .Call('_updogAlpha_obj_offspring', PACKAGE = 'updogAlpha', ocounts, osize, ploidy, prob_geno, bias_val, seq_error, od_param, outlier, out_prop, out_mean, out_disp)
}

#' Just a reparameterization of \code{\link{obj_offspring}}.
#'
#' @inheritParams obj_offspring_vec
#' @param s Same as \code{log(bias_val)} in \code{\link{obj_offspring}}.
#' @param ell We have \code{seq_error = expit(ell)} from \code{\link{obj_offspring}}.
#' @param r Same as \code{log((1.0 - od_param) / od_param)} from \code{\link{obj_offspring}}.
#'
#' @author David Gerard
#'
obj_offspring_reparam <- function(ocounts, osize, ploidy, prob_geno, s, ell, r, outlier = FALSE, out_prop = 0.01) {
    .Call('_updogAlpha_obj_offspring_reparam', PACKAGE = 'updogAlpha', ocounts, osize, ploidy, prob_geno, s, ell, r, outlier, out_prop)
}

#' Same thing as \code{\link{obj_offspring}}, but each sample's log-density has a weight.
#'
#' This is mostly used in the EM algorithm.
#'
#' @inheritParams obj_offspring_vec
#' @param weight_vec A vector of numerics between 0 and 1. They don't have to sum to 1.
#'
#' @author David Gerard
#'
#' @export
#'
obj_offspring_weights <- function(ocounts, osize, weight_vec, ploidy, prob_geno, bias_val = 1, seq_error = 0, od_param = 0, outlier = FALSE, out_prop = 0.01, out_mean = 0.5, out_disp = 1.0 / 3.0) {
    .Call('_updogAlpha_obj_offspring_weights', PACKAGE = 'updogAlpha', ocounts, osize, weight_vec, ploidy, prob_geno, bias_val, seq_error, od_param, outlier, out_prop, out_mean, out_disp)
}

#' Reparameterization of \code{\link{obj_offspring_weights}}.
#'
#' It doesn't make any sense to have outlier = true since the weights are just for the EM already.
#'
#' @inheritParams obj_offspring_weights
#' @inheritParams obj_offspring_reparam
#'
#' @author David Gerard
#'
obj_offspring_weights_reparam <- function(ocounts, osize, weight_vec, ploidy, prob_geno, s, ell, r) {
    .Call('_updogAlpha_obj_offspring_weights_reparam', PACKAGE = 'updogAlpha', ocounts, osize, weight_vec, ploidy, prob_geno, s, ell, r)
}

#' Objective function of outlier part in EM step
#'
#' @param ocounts A vector of reads counts of reference allele.
#' @param osize A vector of total reads counts.
#' @param weight_vec The current probability of a point being an outlier.
#' @param out_mean The current mean of the beta.
#' @param out_disp The current overdispersion parameter of the beta.
#'
#' @author David Gerard
#'
outlier_obj <- function(ocounts, osize, weight_vec, out_mean, out_disp) {
    .Call('_updogAlpha_outlier_obj', PACKAGE = 'updogAlpha', ocounts, osize, weight_vec, out_mean, out_disp)
}

#' Derivative of overdispersion parameter with mean already calculated
#'
#' @inheritParams dbeta_dprop
#'
#' @author David Gerard
#'
#' @seealso dbeta_dtau
#'
dbeta_dtau_withxi <- function(x, n, xi, tau) {
    .Call('_updogAlpha_dbeta_dtau_withxi', PACKAGE = 'updogAlpha', x, n, xi, tau)
}

#' The gradient of \code{\link{outlier_obj}}.
#'
#' @inheritParams outlier_obj
#'
#' @author David Gerard
#'
outlier_grad <- function(ocounts, osize, weight_vec, out_mean, out_disp) {
    .Call('_updogAlpha_outlier_grad', PACKAGE = 'updogAlpha', ocounts, osize, weight_vec, out_mean, out_disp)
}

#' Parent contribution to objective function.
#'
#' For each parental read counts you have, add this to
#' \code{{\link{obj_offspring}}} to get the objective function.
#'
#' @inheritParams obj_offspring_vec
#' @param pcounts The number of counts of the reference allele we observe in
#'     a parent.
#' @param psize The total number of reads we observe in a parent.
#' @param pgeno The genotype of the parent.
#' @param weight A double that multiplies to the log-likelihood (objective).
#'     Should only be used when \code{outlier = FALSE} during the M
#'     step of the EM algorithm.
#'
#' @author David Gerard
#'
obj_parent <- function(pcounts, psize, ploidy, pgeno, bias_val = 1, seq_error = 0, od_param = 0, outlier = FALSE, out_prop = 0.01, out_mean = 0.5, out_disp = 1.0 / 3.0, weight = 1.0) {
    .Call('_updogAlpha_obj_parent', PACKAGE = 'updogAlpha', pcounts, psize, ploidy, pgeno, bias_val, seq_error, od_param, outlier, out_prop, out_mean, out_disp, weight)
}

#' A reparameterization of \code{\link{obj_parent}}
#'
#' @inheritParams obj_offspring_reparam
#' @inheritParams obj_offspring_vec
#' @inheritParams obj_parent
#' @param weight A double between 0 and 1. This should only not be 1 when
#'     doing the em algorithm and the weight is the probabiility that
#'     a point being ok.
#'
#' @author David Gerard
#'
obj_parent_reparam <- function(pcounts, psize, ploidy, pgeno, s, ell, r, weight = 1.0, outlier = FALSE, out_prop = 0.01) {
    .Call('_updogAlpha_obj_parent_reparam', PACKAGE = 'updogAlpha', pcounts, psize, ploidy, pgeno, s, ell, r, weight, outlier, out_prop)
}

#' Gradient of \code{\link{obj_parent}} (when \code{outlier = FALSE}) with respect to
#' \code{bias_val}, \code{seq_error}, and \code{od_param}. Basically, this just calculates
#' derivatives of the log beta-binomial density.
#'
#' @inheritParams obj_parent
#' @inheritParams obj_offspring_vec
#'
#' @return A \code{NumericVector} of length three that contains the partial derivatives of
#'     \code{bias_val}, \code{seq_error}, and \code{od_param}, in that order.
#'
#' @author David Gerard
#'
grad_parent <- function(pcounts, psize, ploidy, pgeno, bias_val = 1, seq_error = 0, od_param = 0, weight = 1.0) {
    .Call('_updogAlpha_grad_parent', PACKAGE = 'updogAlpha', pcounts, psize, ploidy, pgeno, bias_val, seq_error, od_param, weight)
}

#' The gradient of \code{\link{obj_parent_reparam}} with respect to \code{s}, \code{ell},
#' and \code{r}, which recall are reparameterizations of the bias, seqeuncing error, and
#' overdispersion, respectively.
#'
#' @inheritParams obj_offspring_reparam
#' @inheritParams obj_offspring_vec
#' @inheritParams obj_parent
#'
#' @return A \code{NumericVector} of length three with the partial derivatives of \code{s},
#'     \code{ell}, and \code{r} in that order.
#'
#' @author David Gerard
#'
grad_parent_reparam <- function(pcounts, psize, ploidy, pgeno, s, ell, r, weight = 1.0) {
    .Call('_updogAlpha_grad_parent_reparam', PACKAGE = 'updogAlpha', pcounts, psize, ploidy, pgeno, s, ell, r, weight)
}

#' Gradient function for \code{\link{dbetabinom_mu_rho_cpp_double}}.
#'
#' @inheritParams obj_parent
#' @param out_mean The mean of the BB.
#' @param out_disp The overdispersion parameter of the BB.
#' @param weight The probability a point is an outlier.
#'
#' @return A vector of length two. The first element of which is the partial derivative
#'     of the outlier mean. The second element of which is the partial
#'     derivative of the outlier overdispersion parameter.
#'
#' @author David Gerard
#'
out_grad_parent <- function(pcounts, psize, out_mean, out_disp, weight) {
    .Call('_updogAlpha_out_grad_parent', PACKAGE = 'updogAlpha', pcounts, psize, out_mean, out_disp, weight)
}

#' E-step for the parents.
#'
#' @inheritParams obj_parent
#' @inheritParams obj_offspring_vec
#' @param d Same as \code{bias_val} in \code{\link{obj_parent}}.
#' @param eps Same as \code{seq_error} in \code{\link{obj_parent}}.
#' @param tau Same as \code{od_param} in \code{\link{obj_parent}}.
#'
#' @author David Gerard
#'
#'
get_parent_outprop <- function(pcounts, psize, ploidy, pgeno, d, eps, tau, out_prop, out_mean, out_disp) {
    .Call('_updogAlpha_get_parent_outprop', PACKAGE = 'updogAlpha', pcounts, psize, ploidy, pgeno, d, eps, tau, out_prop, out_mean, out_disp)
}

#' Returns the probability of seeing the reference allele after including
#' the mapping-bias and the sequencing-error.
#'
#' @param prob A numeric vector. Each element is the proportion of genomes that contain
#'     the reference allele. This should take on
#'     one of the values 0/K, 1/K, ... , K/K, where K is the ploidy of the individual.
#' @param bias The bias parameter. Should be greater than or equal to zero,
#'     though is typically less than 1. A 1 indicates no bias. A value less than one
#'     indicates a bias towards the reference allele. A value greater than 1 indicates
#'     a bias towards the non-reference allele.
#' @param seq_error The sequencing error rate. This should be between 0 and 1.
#'
#' @author David Gerard
#'
pbias <- function(prob, bias, seq_error) {
    .Call('_updogAlpha_pbias', PACKAGE = 'updogAlpha', prob, bias, seq_error)
}

#' A double version of \code{\link{pbias}}.
#'
#' @inheritParams pbias
#'
#' @seealso \code{\link{pbias}}.
#'
#' @author David Gerard
pbias_double <- function(prob, bias, seq_error) {
    .Call('_updogAlpha_pbias_double', PACKAGE = 'updogAlpha', prob, bias, seq_error)
}

#' An Rcpp version of the \code{\link{dbetabinom}}.
#'
#' @param x A numeric vector of the observed counts.
#' @param size A numeric vector of the number of trials.
#' @param alpha_shape A numeric scalar of the alpha parameter underlying beta.
#' @param beta_shape A numeric scalar of the beta parameter of the underlying beta.
#' @param return_log A logical scalar. Should we return the log of the density (\code{TRUE})
#'     or not (\code{FALSE})?
#'
#' @export
#'
#' @author David Gerard
#'
dbetabinom_cpp <- function(x, size, alpha_shape, beta_shape, return_log = FALSE) {
    .Call('_updogAlpha_dbetabinom_cpp', PACKAGE = 'updogAlpha', x, size, alpha_shape, beta_shape, return_log)
}

#' An Rcpp version of \code{\link{dbetabinom_mu_rho}}.
#'
#' @inheritParams dbetabinom_cpp
#' @param mu The mean of the underlying beta.
#' @param rho The overdispersion parameter of the underlying beta.
#'
#' @author David Gerard
#'
dbetabinom_mu_rho_cpp <- function(x, size, mu, rho, return_log = FALSE) {
    .Call('_updogAlpha_dbetabinom_mu_rho_cpp', PACKAGE = 'updogAlpha', x, size, mu, rho, return_log)
}

#' double version of betabinomial density.
#'
#' @param x The number of successes
#' @param size The number of draws
#' @param mu The mean proportion.
#' @param rho The overdispersion parameter.
#' @param return_log A logical. Should we return the log density
#'     (\code{TRUE}) or not (\code{FALSE})?
#'
#' @author David Gerard
#'
dbetabinom_mu_rho_cpp_double <- function(x, size, mu, rho, return_log) {
    .Call('_updogAlpha_dbetabinom_mu_rho_cpp_double', PACKAGE = 'updogAlpha', x, size, mu, rho, return_log)
}

#' Vectorized version of \code{\link[stats]{dhyper}} for C++ implementation.
#'
#' @param x Vector of quantiles representing the number of white balls drawn without replacement from an urn which contains both black and white balls.
#' @param m The number of white balls in the urn.
#' @param n The number of black balls in the urn.
#' @param k The number of balls drawn from the urn.
#'
#' @author David Gerard
#'
dhyper_cpp <- function(x, m, n, k) {
    .Call('_updogAlpha_dhyper_cpp', PACKAGE = 'updogAlpha', x, m, n, k)
}

#' Rcpp function to get genotype probabilities assuming either
#' F1 population given parental genotypes, the allele frequencies assuming
#' Hardy-Weinberg equilibrium, or just a uniform distribution.
#'
#' @param ploidy The ploidy of the species.
#' @param model Do we assume the genotypes are distributed from an F1 population (\code{"f1"}),
#'     an S1 population (\code{"s1"}),
#'     according to Hardy-Weinberg (\code{hw}), or uniformly (\code{"uniform"})?
#' @param p1geno The first parental genotype if \code{model = "f1"}.
#' @param p2geno The second parental genotype if \code{model = "f1"}.
#' @param allele_freq The allele-frequency if \code{model = "hw"}.
#'
#' @author David Gerard
#' @export
#'
get_prob_geno <- function(ploidy, model, p1geno, p2geno, allele_freq) {
    .Call('_updogAlpha_get_prob_geno', PACKAGE = 'updogAlpha', ploidy, model, p1geno, p2geno, allele_freq)
}

#' Rcpp implementation of \code{\link{get_q_array}}.
#'
#' This function will return the segregation proabilities
#'
#' @inheritParams get_q_array
#'
#' @seealso \code{\link{get_q_array}}.
#'
#' @author David Gerard
#'
#' @export
#'
get_q_array_cpp <- function(ploidy) {
    .Call('_updogAlpha_get_q_array_cpp', PACKAGE = 'updogAlpha', ploidy)
}

#' Log-sum-exponential trick that I use all the time.
#'
#' @param xx A matrix whose rows are to be log-sum-exponentiated.
#'
#' @author David Gerard
#'
logsumexp <- function(xx) {
    .Call('_updogAlpha_logsumexp', PACKAGE = 'updogAlpha', xx)
}

#' Gets all possible binomial probabilities for a given ploidy, bias term, and sequencing
#' error rate.
#'
#' @param ploidy The ploidy of the species.
#' @param bias_val The bias parameter. Should be greater than 0. A value of 1 means no bias.
#' @param seq_error The sequencing error rate.
#'
#' @author David Gerard
#'
#' @export
#'
get_pvec <- function(ploidy, bias_val, seq_error) {
    .Call('_updogAlpha_get_pvec', PACKAGE = 'updogAlpha', ploidy, bias_val, seq_error)
}

#' Stupid implementation of colSums because I guess not implemented in Rcpp sugar.
#'
#' @param x A NumericMatrix.
#'
#' @author David Gerard
#'
#'
colSums_cpp <- function(x) {
    .Call('_updogAlpha_colSums_cpp', PACKAGE = 'updogAlpha', x)
}

#' The expit function.
#'
#' @param x A double.
#'
#' @author David Gerard
#'
expit <- function(x) {
    .Call('_updogAlpha_expit', PACKAGE = 'updogAlpha', x)
}
dcgerard/updogAlpha documentation built on May 14, 2019, 3:10 a.m.