bb.mle: Maximum likelihood estimate for beta binomial distributions
In iZID: Identify Zero-Inflated Distributions

Description Usage Arguments Details Value Reference Examples

View source: R/bb.mle.r

calculate maximum likelihood estimate and the corresponding log likelihood value for beta binomial, beta negative binomial, negative binomial and Poisson distributions.

bb.mle(x, n, alpha1, alpha2, lowerbound = 0.01, upperbound = 10000)

bnb.mle(x, r, alpha1, alpha2, lowerbound = 0.01, upperbound = 10000)

nb.mle(x, r, p, lowerbound = 0.01, upperbound = 10000)

poisson.mle(x)

`x`	A vector of count data. Should be non-negative integers.
`n`	An initial value of the number of trials. Must be a positive number, but not required to be an integer.
`alpha1`	An initial value for the first shape parameter of beta distribution. Should be a positive number.
`alpha2`	An initial value for the second shape parameter of beta distribution. Should be a positive number.
`lowerbound`	A lower searching bound used in the optimization of likelihood function. Should be a small positive number. The default is 1e-2.
`upperbound`	An upper searching bound used in the optimization of likelihood function. Should be a large positive number. The default is 1e4.
`r`	An initial value of the number of success before which m failures are observed, where m is the element of x. Must be a positive number, but not required to be an integer.
`p`	An initial value of the probability of success, should be a positive value within (0,1).

bb.mle, bnb.mle, nb.mle and poisson.mle calculate the maximum likelihood estimate of beta binomial, beta negative binomial, negative binomial and Poisson distributions, respectively.

Please NOTE that the arguments in the four functions are NOT CHECKED AT ALL! The user must be aware of their inputs to avoid getting suspicious results.

Suppose that X is a random count variable that only takes non-negative values. If p has a prior distribution beta(alpha1,alpha2) and X follows a binomial distribution b(n,p), then X follows the beta binomial distribution with

P(X=k)=C(n,k)Beta(k+alpha1,n-k+alpha2)/Beta(alpha1,alpha2),

where C(,) is the combination function, Beta(,) is the beta function and beta(,) stands for the beta distribution.

If X stands for the number of failures observed before the rth success, the probability of X taking the value k under the negative binomial distribution equals

P(X=k)=C(k+r-1,k)p^r(1-p)^k,

As in beta binomial distribution, assume the prior distribution of p is beta(alpha1,alpha2). X follows a beta negative binomial distribution if X follows a negative binomial distribution with parameters r and p. The probability density function of a beta negative binomial distribution is defined as:

P(X=k)=Γ(r+k)Beta(r+alpha1,k+alpha2)/Beta(alpha1,alpha2)/Γ(r)/k!,

where Γ represents the Gamma function.

With the only parameter lambda, the probability density function of a Poisson distribution is

P(X=k)=lambda^k exp(-lambda)/k!

The maximum likelihood estimate of all four distributions can be derived by minimizing the corresponding negative log likelihood function. It is easy to deduce the sample estimate of lambda which is equal to the sample mean. However, it is not so straightforward to solve the optimization problems of the other three distributions. Thus, we adopt the optimization algorithm "L-BFGS-B" by calling R basic function optim. Lower and upper bounds on the unknown parameters are required for the algorithm "L-BFGS-B", which are determined by the arguments lowerbound and upperbound. But note that for the estimate of p, the upper bound for searching is essentially 1-lowerbound.

A row vector containing the maximum likelihood estimate of unknown parameters and the corresponding value of log likelihood.

With bb.mle, the following values are returned:

n: the maximum likelihood estimate of n.
alpha1: the maximum likelihood estimate of alpha1.
alpha2: the maximum likelihood estimate of alpha2.
loglik: the value of log likelihood with maximum likelihood estimates plugged-in.

With bnb.mle, the following values are returned:

r: the maximum likelihood estimate of r.
alpha1: the maximum likelihood estimate of alpha1.
alpha2: the maximum likelihood estimate of alpha2.
loglik: the value of log likelihood with maximum likelihood estimates plugged-in.

With nb.mle, the following values are returned:

r: the maximum likelihood estimate of r.
p: the maximum likelihood estimate of p.
loglik: the value of log likelihood with maximum likelihood estimates plugged-in.

With poisson.mle, the following values are returned:

lambda: the maximum likelihood estimate of lambda.
loglik: the value of log likelihood with maximum likelihood estimate plugged-in.

H. Aldirawi, J. Yang, A. A. Metwally (2019). Identifying Appropriate Probabilistic Models for Sparse Discrete Omics Data, accepted for publication in 2019 IEEE EMBS International Conference on Biomedical & Health Informatics (BHI).

x=extraDistr::rbbinom(2000,12,2,4)
bb.mle(x,3,1,1)
x=extraDistr::rbnbinom(2000,8,3,5)
bnb.mle(x, 3.3, 1, 1)
x=stats::rnbinom(2000,size=5,prob=0.3)
nb.mle(x, 7, 0.5)
x=stats::rpois(2000,7)
poisson.mle(x)

         n   alpha1   alpha2   loglik
[1,] 11.99 2.091543 4.170537 17408.19
           r   alpha1   alpha2   loglik
[1,] 128.913 454.7529 576.1732 14217585
            r         p    loglik
[1,] 5.164495 0.3048159 -6351.504
     lambda    loglik
[1,] 6.9305 -4757.755