Problems to obtain the maximum likelihood estimator

Question

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?

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):

library(VGAM)
n = 10
alpha <- 600
beta <- 400

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

hist(y)

Here is the histogram of those observations:

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)))
    })
}

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

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.