R/npnorm.R

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

#' Compute the non-parametric mixing distrbution
#'
#' This function uses the non-parametric method to compute the mixing
#' distribution based on various losses
#' 
#' \code{npnorm.sq} computes the mixing distribution under squared error loss.
#' 
#' \code{npnorm.cvm} computes the mixing distribution under cramer-von Mises loss.
#' 
#' \code{npnorm.ad} computes the mixing distribution under Anderson-Darling loss.
#'
#' @title Distance-based non-parametric method with normal component.
#' @param data the vector of observation to compute the mixing distribution
#' @param stdev the standard deviation of the component density
#' @param mix the initial mixing distribution
#' @param order the digit to round down the observations. If FALSE, no binning is performed.
#' @param maxit the maximum iterations allowed.
#' @param tol the tolerence. Stop if the directional derivative is less than tol
#' @return the nspmix object with ll meaning the loss.
#' @examples
#' data = rnorm(100, c(0, 2))
#' npnorm.sq(data)
#' npnorm.cvm(data)
#' npnorm.ad(data)
#' datalarge = rnorm(1e5, c(0, 2))
#' npnorm.cvm(datalarge, order = -2)
#' npnorm.ad(datalarge, order = -2)
#' @name npnorm
NULL

#' Compute the non-parametric mixing distrbution with components fixed
#'
#' This function uses the non-parametric method to compute the mixing
#' distribution based on various losses with components fixed when the density is normal.
#' 
#' \code{npnormfc.sq} computes the mixing distribution under squared error loss.
#' 
#' \code{npnormfc.cvm} computes the mixing distribution under cramer-von Mises loss
#' 
#' \code{npnormfc.ll} computes the mixing distribution under KL divergence (Maximum likelihood)
#' 
#' \code{npnormfc.ad} computes the mixing distribution under Anderson-Darling loss.
#'
#' @title Mon-parametric method and components fixed when the component is normal.
#' @param data the vector of observation to compute the mixing distribution
#' @param mu0 the vector of fixed support points
#' @param pi0 the vector of probabilities corresponding to support points
#' @param stdev the standard deviation of the component density
#' @param mix the initial mixing distribution
#' @param order the digit to round down the observations. If FALSE, no binning is performed.
#' @param maxit the maximum iterations allowed.
#' @param tol the tolerence. Stop if the directional derivative is less than tol
#' @return the nspmix object with ll meaning the loss.
#' @examples
#' data = rnorm(100, c(0, 2))
#' npnormfc.ll(data, mu0 = 0, pi0 = 0.2)
#' npnormfc.sq(data, mu0 = 0, pi0 = 0.2)
#' npnormfc.cvm(data, mu0 = 0, pi0 = 0.2)
#' npnormfc.ad(data, mu0 = 0, pi0 = 0.2)
#' datalarge = rnorm(1e5, c(0, 2))
#' npnormfc.cvm(datalarge, order = -2)
#' npnormfc.ad(datalarge, order = -2)
#' npnormfc.ll(datalarge, order = -2)
#' @name npnormfc
NULL

#' Estimate the proportion of a support point
#' 
#' Estimate the proportion of a support point
#' 
#' \code{findnpnormfc.ll} uses \code{npnormfc.ll} to compute the target pi0 and the mixing
#' distribution under the KL divergence (maximum likelihood).
#' 
#' \code{findnpnormfc.cvm} uses \code{npnormfc.cvm} to compute the target pi0 and the mixing
#' distribution under the Cramer-von Mises distance. 
#' 
#' \code{findnpnormfc.ad} uses \code{npnormfc.ad} to compute the target pi0 and the mixing
#' distribution under the Anderson-Darling distance.
#' 
#' @title Estimate the proportion of a support point
#' @param data the vector of observations
#' @param val In maximumu likelihood methods, this is the threshold value. 
#' For the distance method, this is the p.value.
#' @param stdev the standard deviation of normal component
#' @param order the digit to round down the observations. If FALSE, no binning is performed.
#' @return A object of class npnorm.
#' @examples 
#' set.seed(1)
#' data = rnorm(1000, c(0, 2))
#' findnpnormfc.ll(data)
#' findnpnormfc.cvm(data)
#' findnpnormfc.ad(data)
#' set.seed(1)
#' datalarge = rnorm(1e5, c(0, 2))
#' findnpnormfc.ll(datalarge, order = -2)
#' findnpnormfc.cvm(datalarge, order = -2)
#' findnpnormfc.ad(datalarge, order = -2)
#' @name findnpnormfc
NULL

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

FDRnpnormfc = function(neg, pos, result){
  # There might be multiple support points at 0 with no reason. Do sum instead.
  (pnorm(neg, sd = result$beta, lower.tail = TRUE) + pnorm(pos, sd = result$beta, lower.tail = FALSE)) * sum(result$mix$pr[result$mix$pt == 0])/
    (pnpnorm1(neg, result$mix$pt, result$mix$pr, stdev = result$beta, lt = TRUE) + pnpnorm1(pos, result$mix$pt, result$mix$pr, stdev = result$beta, lt = FALSE))
}

#' Find the rejection region with normal density
#' 
#' Find the rejection region assuming normal component from the npnormfc 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 normal density
#' @param result an object of npnormfc
#' @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.npnormfc = function(result, alpha = 0.05){
  # object check done via npnormfc2npnorm
  ans = nloptr(c(-1, 1), function(vec, result, alpha){
    v = tan(vec)
    -pnpnorm1(v[1], result$mix$pt, result$mix$pr, result$beta, TRUE) - pnpnorm1(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)
    FDRnpnormfc(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)
}

## MAXIMUM LIKELIHOOD

#' @rdname npnormfc
#' @export
npnormfc.ll = function(data, mu0 = 0, pi0 = 0, stdev = NULL, mix = NULL, order = FALSE, maxit = 1000, tol = 1e-6){
  if (is.null(stdev)) stdev = 1
  if (is.null(mix)) mix = initnpnorm(data, mu0, pi0, stdev)
  if (!order){
      r = npnormfcll(data = data, mu0fixed = mu0, pi0fixed = pi0, stdev = stdev,
             mu0_ = mix$pt, pi0_ = mix$pr, maxit = maxit, tol = tol)
    }else{
      x = bin(data, order = order)
      r = npnormfcllweighted(data = x$v, weights = x$w, mu0fixed = mu0, pi0fixed = pi0, 
        stdev = stdev, h = 10^order, mu0_ = mix$pt, pi0_ = mix$pr, maxit = maxit, tol = tol)
    }
    attr(r, "class") = "nspmix" 
    r
}

#' @rdname findnpnormfc
#' @export
findnpnormfc.ll = function(data, val = log(length(data)), stdev = 1, order = FALSE){
  if (!order){
    r = findnpnormll(data, val, stdev)
  }else{
    x = bin(data, order = order)
    r = findnpnormllweighted(data = x$v, weights = x$w, val = val, stdev = stdev, h = 10^order)
  }
  attr(r, "class") = "nspmix" 
  r
}

## SQUARED ERROR

#' @rdname npnorm
#' @export
npnorm.sq = function(data, stdev = NULL, mix = NULL, order = FALSE, maxit = 100, tol = 1e-6){
  if (is.null(stdev)) stdev = 1
  if (is.null(mix)) mix = initnpnorm(data, 0, 0, stdev)
  if (!order){
    r = npnormsq(data, stdev, mix$pt, mix$pr, maxit, tol)     
  }else{
    stop("Not yet implemented")
  }  
  attr(r, "class") = "nspmix" 
  r
}

#' @rdname npnormfc
#' @export
npnormfc.sq = function(data, mu0 = 0, pi0 = 0, stdev = NULL, mix = NULL, order = FALSE, maxit = 100, tol = 1e-6){
  if (is.null(stdev)) stdev = 1
  if (is.null(mix)) mix = initnpnorm(data, mu0, pi0, stdev)
  if (!order){
      r = npnormfcsq(data, mu0, pi0, stdev, mix$pt, mix$pr, maxit, tol)     
  }else{
    stop("Not yet implemented")
  } 
  attr(r, "class") = "nspmix" 
  r
}

## CRAMER-VON MISES

#' @rdname npnorm
#' @export
npnorm.cvm = function(data, stdev = NULL, mix = NULL, order = FALSE, maxit = 100, tol = 1e-6){
  if (is.null(stdev)) stdev = 1
  if (is.null(mix)) mix = initnpnorm(data, 0, 0, stdev)
  if (!order){
    r = npnormcvm(data, stdev, mix$pt, mix$pr, maxit, tol)     
  }else{
    x = bin(data, order = order)
    r = npnormcvmweighted(x$v, stdev, 10^order, mix$pt, mix$pr, x$w, maxit, tol)
  }   
  attr(r, "class") = "nspmix" 
  r
}

#' @rdname npnormfc
#' @export
npnormfc.cvm = function(data, mu0 = 0, pi0 = 0, stdev = NULL, mix = NULL, order = FALSE, maxit = 100, tol = 1e-6){
  if (is.null(stdev)) stdev = 1
  if (is.null(mix)) mix = initnpnorm(data, mu0, pi0, stdev)
  if (!order){
    r = npnormfccvm(data, mu0, pi0, stdev, mix$pt, mix$pr, maxit, tol)     
  }else{
    x = bin(data, order = order)
    r = npnormfccvmweighted(x$v, mu0, pi0, stdev, 10^order, mix$pt, mix$pr, x$w, maxit, tol)
  } 
  attr(r, "class") = "nspmix" 
  r
}

#' @rdname findnpnormfc
#' @export
findnpnormfc.cvm = function(data, val = 0.05, stdev = 1, order = FALSE){
  if (!order){
      r = findnpnormcvm(data, qcvm(val, lower.tail = FALSE), stdev)      
  }else{
      x = bin(data, order = order)
      r = findnpnormcvmweighted(x$v, qcvm(val, lower.tail = FALSE), stdev, 10^order, x$w)
  } 
  attr(r, "class") = "nspmix" 
  r
}

## ANDERSON-DARLING

#' @rdname npnorm
#' @export
npnorm.ad = function(data, stdev = NULL, mix = NULL, order = FALSE, maxit = 100, tol = 1e-6){
  if (is.null(stdev)) stdev = 1
  if (is.null(mix)) mix = initnpnorm(data, 0, 0, stdev)
  if (!order){
    r = npnormad(data, stdev, mix$pt, mix$pr, maxit, tol)  
  }else{
    x = bin(data, order = order)
    r = npnormadweighted(x$v, stdev, 10^order, mix$pt, mix$pr, x$w, maxit, tol)
  }
  attr(r, "class") = "nspmix" 
  r
}

#' @rdname npnormfc
#' @export
npnormfc.ad = function(data, mu0 = 0, pi0 = 0, stdev = NULL, mix = NULL, order = FALSE, maxit = 100, tol = 1e-6){
  if (is.null(stdev)) stdev = 1
  if (is.null(mix)) mix = initnpnorm(data, mu0, pi0, stdev)
  if (!order){
    r = npnormfcad(data = data, mu0fixed = mu0, pi0fixed = pi0, stdev = stdev,
                   mu0_ = mix$pt, pi0_ = mix$pr, maxit = maxit, tol = tol)
  }else{
    x = bin(data, order = order)
    r = npnormfcadweighted(x$v, mu0, pi0, stdev, 10^order, mix$pt, mix$pr, x$w, maxit, tol)
  }
  attr(r, "class") = "nspmix" 
  r
}

#' @rdname findnpnormfc
#' @export
findnpnormfc.ad = function(data, val = 0.05, stdev = 1, order = FALSE){
  if (!order){
    r = findnpnormad(data, qad(val, lower.tail = FALSE), stdev)  
  }else{
    x = bin(data, order = order)
    r = findnpnormadweighted(x$v, qad(val, lower.tail = FALSE), stdev, 10^order, x$w)
  }
  attr(r, "class") = "nspmix" 
  r
}