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)
A vector of count data. Should be non-negative integers.
An initial value of the number of trials. Must be a positive number, but not required to be an integer.
An initial value for the first shape parameter of beta distribution. Should be a positive number.
An initial value for the second shape parameter of beta distribution. Should be a positive number.
A lower searching bound used in the optimization of likelihood function. Should be a small positive number. The default is 1e-2.
An upper searching bound used in the optimization of likelihood function. Should be a large positive number. The default is 1e4.
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.
An initial value of the probability of success, should be a positive value within (0,1).
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).
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 \(r\)th 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)=\Gamma(r+k)Beta(r+alpha1,k+alpha2)/Beta(alpha1,alpha2)/\Gamma(r)/k!\),
where \(\Gamma\) 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.
# NOT RUN {
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)
# }
Run the code above in your browser using DataLab