R/sample.h.r

In iZID: Identify Zero-Inflated Distributions

#' Generate random deviates from zero-inflated or hurdle models
#'
#' Generate random deviates from zero-inflated or hurdle Poisson, negative binomial, beta binomial and beta negative binomial
#' models.
#'
#' By setting \code{distri=poisson}, \code{sample.h} and \code{sample.zi} simulates \code{N} random deviates from hurdle and
#' zero-inflated Poisson distribution, respectively, and so on forth. For arguments with length larger than 1, only the first
#' element will be used.
#'
#' Arguments \code{r} and \code{p} are for the use of zero-inflated and hurdle negative binomial distributions. \code{alpha1},
#' \code{alpha2} and \code{n} are for
#' zero-inflated and hurdle beta binomial distributions. \code{r}, \code{alpha1} and \code{alpha2} are used in zero-inflated and hurdle beta
#' negative binomial distributions.
#'
#' The procedure of generating random deviates follows the work of Aldirawi et al. (2019). The algorithm calls functions for
#' standard distributions to simulate the non-zero samples. Random deviates from standard Poisson and negative binomial
#' distributions can be generated by basic R function \code{rpois} and \code{rnbinom}. Functions rbbinom and rbnbinom are
#' available for standard beta binomial and beta negative binomial distributions in R package \code{extraDistr}.
#'
#' @param N  The sample size. Should be a positive number. If it is not an integer, N will be automatically rounded up to the
#' smallest integer that no less than N.
#' @param phi The structural parameter \eqn{\phi}, should be a positive value within (0,1).
#' @param distri The corresponding standard distribution. Can be one of \{'poisson','nb','bb', 'bnb'\}, which corresponds to Poisson,
#' negative binomial, beta binomial and beta negative binomial distributions respectively.
#' @param lambda A value for the parameter of Poisson distribution. Should be a positive number.
#' @param r the number of success before which m failures are observed, where m is a random variable from negative binomial
#' or beta negative binomial distribution. Must be a positive number. If it is not an integer, r will be automatically rounded up to the
#' smallest integer that no less than r.
#' @param p The probability of success, should be a positive value within (0,1).
#' @param alpha1 The first shape parameter of beta distribution. Should be a positive number.
#' @param alpha2 The second shape parameter of beta distribution. Should be a positive number.
#' @param n The number of trials. Must be a positive number. If it is not an integer, n will be automatically rounded up to the
#' smallest integer that no less than n.
#'
#' @return A vector of length \eqn{N} containing non-negative integers from the zero-inflated or hurdle version of distribution
#' determined by \eqn{distri}.
#' @export
#'
#'@section Reference:
#'\itemize{
#'  \item{H. Aldirawi, J. Yang, A. A. Metwally, Identifying Appropriate Probabilistic Models for Sparse Discrete Omics Data,
#'  accepted for publication in 2019 IEEE EMBS International Conference on Biomedical & Health Informatics (BHI) (2019).}
#'  \item{T. Wolodzko, extraDistr: Additional Univariate and Multivariate Distributions, R package version 1.8.11 (2019),
#'   https://CRAN.R-project.org/package=extraDistr.}
#'}
#'
#' @examples
#' t1=sample.h(N=2000,phi=0.2,distri='Poisson',lambda=5)  ##hurdle poisson random values
#' t2=sample.h(N=2000,phi=0.2,distri='nb',r=10,p=0.6)   ##hurdle negative binomial
#' t3=sample.h(N=2000,phi=0.2,distri='bb',alpha1=8,alpha2=9,n=10)   ##hurdle beta binomial
#' ##hurdle beta negative binomial.
#' t4=sample.h(N=2000,phi=0.2,distri='bnb',r=10,alpha1=8,alpha2=9)
#'
#'t1=sample.zi(N=2000,phi=0.2,distri='Poisson',lambda=5)  ##zero-inflated poisson random values
#' t2=sample.zi(N=2000,phi=0.2,distri='nb',r=10,p=0.6)   ##zero-inflated negative binomial
#' t3=sample.zi(N=2000,phi=0.2,distri='bb',alpha1=8,alpha2=9,n=10)   ##zero-inflated beta binomial
#' ##zero-inflated beta negative binomial
#' t4=sample.zi(N=2000,phi=0.2,distri='bnb',r=10,alpha1=8,alpha2=9)
sample.h<-function(N,phi,distri='poisson',lambda=NA,r=NA,p=NA,alpha1=NA,alpha2=NA,n=NA){
  ##check the validity of some inputs
  distri=tolower(distri)[1]  ##so that all letters are in lower case. only use the first value if user inputs a vector but not a scalar.
  N=ceiling(N)[1]
  phi=phi[1]
  if(N<=0)
    stop('the sample size N is too small.')
  if((phi>=1)|(phi<=0))
    stop('phi is not between 0 and 1 (not including 0 and 1).')
  if(!(distri%in%c('poisson','nb','bb','bnb')))
   stop('please input a distribution name among poisson,nb,bb,bnb.')
  if(distri=="poisson"){
    lambda=lambda[1]
	if(lambda<=0)
	   stop('lambda is less or equal to 0.')
  }
  if(distri=="nb"){
    r=ceiling(r[1])
	p=p[1]
	if(r<=0)
	   stop('r is too small.')
	if((p<=0)|(p>=1))
	   stop('p is not between 0 and 1 (not including 0 and 1).')
  }
  if((distri=="bb")|(distri=="bnb")){
    alpha1=alpha1[1]
	alpha2=alpha2[1]
	if(alpha1<=0)
	   stop('alpha1 is less or equal to 0.')
	if(alpha2<=0)
	   stop('alpha2 is less or equal to 0.')
  }
  if(distri=="bb"){
    n=ceiling(n[1])
	if(n<=0)
	   stop('n is too small.')
  }
  if(distri=="bnb"){
    r=ceiling(r[1])
	if(r<=0)
	   stop('r is too small.')
  }

  ##generate random values from hurdle models
  ans=stats::rbinom(N,size=1,prob=1-phi)
  m=sum(ans==1)
  p0=switch(distri,nb=p^r,bb=beta(alpha1,n+alpha2)/beta(alpha1,alpha2),bnb=beta(alpha1+r,alpha2)/beta(alpha1,alpha2),exp(-lambda)) ##prob of z=0
  M=ceiling((m+2*sqrt(p0*(m+p0))+2*p0)/(1-p0))
  z=switch(distri,nb=stats::rnbinom(M,size=r,prob=p),bb=extraDistr::rbbinom(M,n,alpha1,alpha2),bnb=extraDistr::rbnbinom(M,r,alpha1,alpha2),stats::rpois(M,lambda))
  u=z[z>0]
  t=length(u)
  if(t<m){
    u1=NULL
    itemp=m-t
    while(itemp>0){
      temp=switch(distri,nb=stats::rnbinom(itemp,size=r,prob=p),bb=extraDistr::rbbinom(itemp,n,alpha1,alpha2),bnb=extraDistr::rbnbinom(itemp,r,alpha1,alpha2),stats::rpois(itemp,lambda))
      u1=c(u1,temp[temp>0])
	  itemp=m-t-length(u1)
    }
    u=c(u,u1)
  }
  ans[ans==1]=u[1:m]
  return(ans)
}