R/nptfixedcomp.R

In xiangjiexue/npfixedcomp: Non-parametric methods for estimating mixing distribution

#' load npnormfc functions
#' 
#' load npnormfc functions
#' 
#' @title load npnormfc functions
loadnptfc = function(){
  source(system.file("extdata", "nptfc.R", package = "npfixedcomp"))
}

#' Find the rejection region with t density
#' 
#' Find the rejection region assuming t component from the nptfc object.
#' The rejection region is calculated using the density estimate rather than data points hence robust.
#' The rejection is based on the hypothesis is located at 0.
#' The optimisation is done via NLopt library (The package nloptr)
#' 
#' @title Find the rejection region with t density
#' @param result an object of nptfc
#' @param alpha the FDR controlling rate.
#' @return a list with par is the boundary for rejection and area is the propotion of rejection
#' @export
rejectregion.nptfc = function(result, alpha = 0.05){
  # object check done via nptfc2npt
  ans = nloptr(c(-1, 1), function(vec, result, alpha){
    v = tan(vec)
    -pnpt1(v[1], result$mix$pt, result$mix$pr, result$beta, TRUE) - pnpt1(v[2], result$mix$pt, result$mix$pr, result$beta, FALSE)
  }, lb = c(-base::pi / 2, 0), ub = c(0, base::pi / 2), eval_g_ineq = function(vec, result, alpha){
    v = tan(vec)
    FDRnptfc(v[1], v[2], result) - alpha
  }, result = result, alpha = alpha, opts = list(algorithm = "NLOPT_GN_ORIG_DIRECT", maxeval = 10000))
  list(par = tan(ans$solution), area = -ans$objective)
}

FDRnptfc = function(neg, pos, result){
  (pt(neg, df = result$beta, lower.tail = TRUE) + pt(pos, df = result$beta, lower.tail = FALSE)) * sum(result$mix$pr[result$mix$pt == 0])/
    (pnpt1(neg, result$mix$pt, result$mix$pr, stdev = result$beta, lt = TRUE) + pnpt1(pos, result$mix$pt, result$mix$pr, stdev = result$beta, lt = FALSE))
}

#' compute the NPMLE with t component with fixed components and probabilities
#' 
#' compute the NPMLE with t component with fixed components and probabilities
#' 
#' @title NPMLE with t component with fixed components and probabilities
#' @param data the vector of observations to compute NPMLE
#' @param mu0 the vector of fixed support points
#' @param pi0 the vector of probabilities corresponding to support points
#' @param df the degree of freedom of component t distribution
#' @param mix the initial mixing distribution
#' @param maxit the maximum iteration allowed.
#' @param tol the tolarence. 
#' @return a object of class `npt`
#' @rdname nptfcll
#' @export
nptfc.ll = function(data, mu0 = 0, pi0 = 0, df, mix = NULL, maxit = 100, tol = 1e-6){
  if (is.null(mix)) mix = initnpnorm(data, mu0, pi0, sqrt(df / (df - 2)))
  nptfcll(data = data, mu0fixed = mu0, pi0fixed = pi0, stdev = df,
          mu0_ = mix$pt, pi0_ = mix$pr, maxit = maxit, tol = tol)
}

#' Estimate the proportion of 0 using likelihood-based method for t component
#' 
#' Estimate the proportion of 0 using likelihood-based method for t component
#' 
#' 
#' @title Estimate the proportion of 0 using likelihood-based method
#' @param data the vector of observations
#' @param val the threshold
#' @param df the degree of freedom
#' @return A list with result an object of family nptfc with estimated pi0 and optim.stat
#' containing the information using optim. The full mixing distribution can be
#' found using nptfc2npt
#' @rdname findnptfcll
#' @export
findnptfc.ll = function(data, val = log(length(data)), df){
  findnptll(data = data, val = val, stdev = df)  
}

#' density of t distribution
#' 
#' Computes the t distribution with primary focus on evaluating non-central t distribution
#' 
#' Two methods are provided for evaluating the density. The first kind is using the numerical
#' intergtation. The result should be precise, but it uses more time when df is low. Use this if
#' you want demand a high accuracy. The implementation for |ncp| < 5 and |x| < 5 and/or 
#' |ncp-x|<10.0 is provided by R::dnt.
#' 
#' The second one is of the integration by part type. This is fast and relative accurate when
#' df is small. when df is large use numerical integration or normal integration. Using this 
#' method, the density is evaluated using sucessive number of gaussian integral.
#' 
#' @title Evaluating density of t distribution
#' @param x the vector of quantiles.
#' @param df the degree of freedom. can be real if use.integration = TRUE, integer otherwise.
#' @param ncp the non-centrality parameters.
#' @param log logical; if TRUE, the density is return on log scale.
#' @param use.integration logical; if TRUE, the density is evaluated using numerical integration.,
#' @return a matrix of length(x) * length(ncp)
#' @export
dt = function(x, df, ncp, log = FALSE, use.integration = FALSE){
  if (df > 20 & use.integration){
    warning("df is too big for default method, return results using integration.")
    use.integration = TRUE
  }
  dt1(x, df, ncp, log, use.integration)
}