I'm trying to calcul the MLE from a Beta-Binomial distribution. However, I'm having problems defining the estimator.

I'm using the following function:

$$ \widehat{\ell\,}(\theta\,;x)=\sum_{i=1}^n \ln f(x_i\mid\theta), $$ I tried to calcul the function and then apply the log() function in each point:

n = 9
k = seq(0,n)

y = VGAM::dbetabinom(x = k, size = n, prob = 0.09894759, rho = 0.08603263)

mle = sum(log(y))

The problem is that the results I get are different from when compared to some function already existing in R.

I'm using the function below and trying to estimate alpha e beta from it. The parameters n and k are predefined.

$$ f(k\mid n,\alpha,\beta) = \frac{\Gamma(n+1)}{\Gamma(k+1)\Gamma(n-k+1)} \frac{\Gamma(k+\alpha)\Gamma(n-k+\beta)}{\Gamma(n+\alpha+\beta)} \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} $$

I imagine that I have a problem with the concept, but I can not find the correct result. I don't undersand why when I sum all log() values I don't get the MLE.

How do I make the correct calculation of this estimator?

  • $\begingroup$ Where do you optimise the (log-)likelihood? $\endgroup$
    – Xi'an
    Jul 1, 2019 at 15:28
  • $\begingroup$ Actually, I'm trying to optimize this function with Excel's Solver, but when I compare my results to those obtained by functions in R, they are different. $\endgroup$ Jul 2, 2019 at 7:32

1 Answer 1


You are quite confused. I see two major confusions: you are confusing the idea of the log-likelihood with that of the log-density and you do not understand the idea of an estimator.

The density function of a beta binomial is given by the formula you listed and can be evaluated in R by the function dbetabinom although note that dbetabinom actually uses a different parameterization than the formula you listed. You can get the same parameterization if you use dbetabinom.ab.

The log-likelihood is a function of the parameters ($\alpha$, and $\beta$) given some observed data. So suppose you have 1000 independent observations, $k_1, \ldots,k_{1000}$ and you assume each one came from a BetaBinomial with the same $\alpha$ and $\beta$. This is usually written $$ k_i \sim \text{BetaBinomial}(n, \alpha, \beta) $$

In R we can draw such a sample of independent observations from the BetaBinomial like this (sorry I used y in my code rather than k):

n = 10
alpha <- 600
beta <- 400

y <- rbetabinom.ab(1000, size=n, alpha, beta)


Here is the histogram of those observations:

enter image description here

The log-likelihood, $\ell(\alpha, \beta)=\log{L(\alpha, \beta)}$ is the log of the likelihood function which is a function of the parameters. Because we assumed each observation is independent we can factorize this likelihood as

$$ L(\alpha, \beta) = \prod_{i=1}^{1000} p(k_i|n, \alpha, \beta) $$ where $p(k_i|n,\alpha,\beta)$ is the formula you listed.

The MLE (maximum likelihood estimator) is the $\alpha$ and $\beta$ which maximize this log-likelihood function

$$ (\hat{\alpha}, \hat{\beta}) = \displaystyle\text{argmax}_{\alpha, \beta}{\ell(\alpha, \beta)} $$

Now for many distributions you can get a closed form solution for the MLE. There is no closed form solution for the MLE of a beta-binomial. So you would have to obtain the MLE using computational means (Newton-Raphson algorithm etc).

Now for illustrative purposes Lets look at the log-likelihood function with one of parameters held constant at its true value:

f1 <- function(a){
  sapply(a, function(ax){ sum(log(dbetabinom.ab(y, size=n, ax, beta)))})

f2 <- function(b){
  sapply(b, function(bx){ 
    sum(log(dbetabinom.ab(y, size=n, alpha, bx)))

curve(f1(x), 0, 1000)
abline(v=alpha, col="gray", lty=2)
curve(f2(x), 0, 1000)
abline(v=beta, col="gray", lty=2)

enter image description here

f1 is $\ell(\alpha, \beta=400)$ and similarly f2 is $\ell(\alpha=600, \beta)$ I have plotted the log-likelihood as these parameters vary and in a vertical gray dotted line I have the true value. You can see this true value is close to the value which would maximize these functions.

  • $\begingroup$ I understood your example. What I'm trying to do is optimize this function using Excel's Solver. For this I need to select the objective cell and the variable cells. I already did the calculation using the RSS estimator (I select the RSS value as the objective and the alpha and beta parameters as variables). However, I can not define a cell with the MLE value to use as an objective cell in Excel. $\endgroup$ Jul 2, 2019 at 7:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.