R/RcppExports.R

In ebTobit: Empirical Bayesian Tobit Matrix Estimation

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

#' Nonparametric Maximum Likelihood via Expectation Maximization
#'
#' Compute the nonparametric maximum likelihood estimate given a likelihood
#' matrix. The matrix A is structured so that A_\{ij\} = f(X_i | theta_j) for
#' some grid of potential parameter values theta_1, ..., theta_p and
#' observations X_1, ..., X_n. The parameters, theta_j, can be multidimensional
#' because all that is required is the likelihood. Convergence is achieved when
#' the relative improvements of the log-likelihood is below the provided
#' tolerance level.
#'
#' @param A numeric matrix likelihoods
#' @param maxiter early stopping condition
#' @param rtol convergence tolerance: abs(loss_new - loss_old)/abs(loss_old)
#' @return the estimated prior distribution (a vector of masses corresponding
#' to the columns of A)
#'
#' @examples
#' set.seed(1)
#' t = sample(c(0,5), size = 100, replace = TRUE)
#' x = t + stats::rnorm(100)
#' gr = seq(from = min(x), to = max(x), length.out = 50)
#' A = stats::dnorm(outer(x, gr, "-"))
#' EM(A)
#' \dontrun{
#' # compare to solution from rmosek (requires additional library installation):
#' all.equal(
#'     REBayes::KWPrimal(A = A, d = rep(1, 50), w = rep(1/100, 100))$f,
#'     EM(A, maxiter = 1e+6, rtol = 1e-16), # EM alg converges slowly
#'     tolerance = 0.01
#' )
#' }
#' @useDynLib ebTobit
#' @importFrom Rcpp evalCpp
#' @export
EM <- function(A, maxiter = 1e+4L, rtol = 1e-6) {
    .Call('_ebTobit_EM', PACKAGE = 'ebTobit', A, maxiter, rtol)
}

#' Helper Function - generate likelihood for pair (L,R) and mean gr
#'
#' Compute P(L_i, R_i | theta = t_k) for observations (L_i, R_i) and grid of
#' mean t_k.
#'
#' @param L numeric vector of lower bounds
#' @param R numeric vector of upper bounds
#' @param gr numeric vector of means
#' @param s1 numeric vector of standard deviations
#' @return the likelihood under partial interval censoring
#'
#' @examples
#' # set-up
#' p = 15
#' gr = stats::rnorm(p)
#' L = R = stats::rnorm(p)
#' missing.idx = sample.int(n = p, size = p/5)
#' L[missing.idx] = L[missing.idx] - stats::runif(length(missing.idx), 0, 1)
#' R[missing.idx] = R[missing.idx] + stats::runif(length(missing.idx), 0, 1)
#'
#' # R solution
#' lik = prod(ifelse(
#'            L == R,
#'            stats::dnorm(L-gr),
#'            stats::pnorm(R-gr) - stats::pnorm(L-gr)))
#'
#' # Compare R to RcppParallel method
#' all.equal(lik, lik_GaussianPIC(L, R, gr, rep(1,p)))
#' @useDynLib ebTobit
#' @importFrom Rcpp evalCpp
#' @export
lik_GaussianPIC <- function(L, R, gr, s1) {
    .Call('_ebTobit_lik_GaussianPIC', PACKAGE = 'ebTobit', L, R, gr, s1)
}

#' Helper Function - generate likelihood matrix
#'
#' Compute a matrix L whose entries are L[i,k] = P(L_i, R_i | theta = t_k) for
#' observations (L_i, R_i) and grid of means t_k.
#'
#' @param L n x p matrix of lower bounds
#' @param R n x p matrix of upper bounds
#' @param gr m x p matrix of candidate means
#' @param s1 n x p matrix of standard deviations
#' @return the n x m likelihood matrix under partial interval censoring
#'
#' @examples
#' # set-up
#' n = 100; m = 50; p = 5
#' gr = matrix(stats::rnorm(m*p), m, p)
#' L = R = matrix(stats::rnorm(n*p), n, p)
#' s1 = matrix(1, n, p)
#' missing.idx = sample.int(n = n*p, size = p*p)
#' L[missing.idx] = L[missing.idx] - stats::runif(p, 0, 1)
#'
#' # R solution
#' lik = matrix(nrow = n, ncol = m)
#' for (i in 1:n) {
#'     for(k in 1:m) {
#'         lik[i,k] = prod(ifelse(
#'             L[i,] == R[i,],
#'             stats::dnorm(L[i,]-gr[k,], sd = s1[i,]),
#'             stats::pnorm(R[i,]-gr[k,], sd = s1[i,]) - stats::pnorm(L[i,]-gr[k,], sd = s1[i,])
#'         ))
#'     }
#' }
#'
#' # Compare R to RcppParallel method
#' all.equal(lik, likMat(L, R, gr, s1))
#' @useDynLib ebTobit
#' @importFrom Rcpp evalCpp
#' @import RcppParallel
#' @export
likMat <- function(L, R, gr, s1) {
    .Call('_ebTobit_likMat', PACKAGE = 'ebTobit', L, R, gr, s1)
}