R/npfixedcomp.R

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

#' npfixedcomp
#' 
#' npfixedcomp
#' 
#' @docType package
#' @author Xiangjie Xue
#' @import Rcpp RcppEigen RcppNumerical
#' @importFrom nloptr nloptr
#' @importFrom stats pnorm pt uniroot rnorm
#' @useDynLib npfixedcomp
#' @name npfixedcomp
NULL

#' Bin the continuous data.
#' 
#' This function bins the continuous data using the equal-width bin in order to 
#' speed up some functions in this package.
#' 
#' h is taken as 10^(-order)
#'
#' the observations are rounded down to the bins ..., -h, 0, h, ...
#' 
#' To further speed up the process, omit the bin that has 0 count.
#' 
#' @title Bin the continuous data.
#' @param data the observation to be binned.
#' @param order see the details
#' @return a list with v be the representative value of each bin and w be the count in the corresponding bin.
#' @export
bin = function(data, order = -5){
  h = 10^order
  data = floor(data / h) * h
  t = table(data)
  index = t != 0
  list(v = as.numeric(names(t))[index], w = as.numeric(t)[index])
}

# The p function of ad distribution.c
pad = function(z, lower.tail = TRUE){
  ans = numeric(length(z))
  ans[z == 0] = 0
  ans[z == Inf] = 1
  index = z != 0 & z != Inf
  j = 0;
  res = choose(-0.5, j) * (4 * j + 1) * exp(-(4 * j + 1)^2 * base::pi^2 / 8 / z[index]) * adintegralvec(z[index], j)
  ans[index]= res;
  while (sum(res^2) > 1e-10){
    j = j + 1; 
    res = choose(-0.5, j) * (4 * j + 1) * exp(-(4 * j + 1)^2 * base::pi^2 / 8 / z[index]) * adintegralvec(z[index], j)
    ans[index] = ans[index] + res;
  }
  ans[index] = ans[index] * sqrt(2 * base::pi) / z
  if (lower.tail){
    ans
  }else{
    1 - ans
  }
}

qad = function(p, lower.tail = TRUE){
  sapply(p, function(ddd){
    uniroot(function(l, lower.tail) {pad(l, lower.tail = lower.tail) - ddd},
            interval = c(0, 500), lower.tail = lower.tail, 
            extendInt = ifelse(lower.tail, "upX", "downX"))$root
  })
}

# The a1(z) in Anderson and Darling (1952). The asymptotic distribution function of cvm statistics
# it was mentioned in the paper that u_{j + 1} / u{j} < 1.0007exp(-(4j + 1)/2/z)
pcvm = function(x, lower.tail = TRUE){
  ans = numeric(length(x))
  ans[x == 0] = 0
  ans[x == Inf] = 1
  index = x != 0 & x != Inf
  # calculate the number of terms
  if (sum(index) > 0){
    jmax = ceiling((log(1e-10) * max(x[index]) * -2 - 1) / 4)
    x.pt = rep(x[index], rep(jmax + 1, length(x[index])))
    j = 0:jmax
    ans[index] = .colSums((-1)^j * choose(-0.5, j) * sqrt(4 * j + 1) * exp(-(4 * j + 1)^2/16/x.pt) * 
                            besselK((4 * j + 1)^2 / 16/ x.pt, nu = 1/4), m = jmax + 1, n = length(x[index])) / 
      base::pi / sqrt(x[index])
  }
  if (lower.tail) {ans}else{1 - ans}
}

qcvm = function(p, lower.tail = TRUE){
  sapply(p, function(ddd){
    uniroot(function(l, lower.tail) {pcvm(l, lower.tail = lower.tail) - ddd},
            interval = c(0, 500), lower.tail = lower.tail, 
            extendInt = ifelse(lower.tail, "upX", "downX"))$root
  })
}

#' plot npnormdist object
#' 
#' plot npnormdist object. Plotting is done via plot.npnorm.
#' mask the ll element since the ll in npnorm is the likelihood
#' while the ll element in npnormdist is the loss under cvm.
#' 
#' @title plot npormdist object
#' @param result an object of family npnormdist
#' @param ... other argument passed to the underlying functions
#' @export
plotnpnormdist = function(result, ...){
  if (result$family != "npnormdist"){
    stop("Input is not the family of npnormdist")
  }
  result$family = "npnorm"
  result$ll = NULL
  plot(result, ...)
}