"How to solve a system of equations of estimators (alpha and beta) using the method of maximum likelihood (beta distribution) in R?"

Question

set.seed(42)
# Data
x <- rbeta(100, 4.0, 2.0)

logLikGrad <- function(params){
 a <- params[1]
 b <- params[2]
 
 grad_a <- (digamma(a+b)) - (digamma(a)) + sum(log(x))
 grad_b <- (digamma(a+b)) - (digamma(b)) + sum(log(1-x))
 
 
 c(-grad_a, -grad_b)
}


startParams <- c(4, 2)


fit <- nleqslv(startParams, logLikGrad)


fit$x

Hello, I'm encountering an issue with my code. It's yielding negative values that are far from correct. I've searched online and even consulted ChatGPT, but haven't found a solution to my problem. Could anyone help? The task was to derive estimators for the parameters α and β using the method of maximum likelihood (ML) and solve it, implementing it in RStudio.

Answer 1

Apart from a one significant issue plus a minor point, your code works.

The minor point is that if you want other people to reproduce your results then you need to specify what packages you are using, here nleqslv

The significant issue is that you started with a likelihood of $$\prod_i \left(\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}x_i^{\alpha-1}(1-x_i)^{\beta-1}\right) = \left(\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\right)^n\left(\prod_i x_i\right)^{\alpha-1}\left(\prod_i (1-x_i)\right)^{\beta-1}$$ but, in finding the log-likelihood and its derivatives, you lost the $n$ term. You need to include it.

So if your code looked something like

library(nleqslv)
set.seed(42)
# Data
x <- rbeta(100, 4.0, 2.0)

logLikGrad <- function(params){
 a <- params[1]
 b <- params[2]
 
 grad_a <- length(x) * (digamma(a+b) - digamma(a)) + sum(log(x))
 grad_b <- length(x) * (digamma(a+b) - digamma(b)) + sum(log(1-x))
 
 
 c(-grad_a, -grad_b)
}

startParams <- c(4, 2)


fit <- nleqslv(startParams, logLikGrad)


fit$x
 

then it would give the results [1] 3.828875 1.992556 which are not far from the original values of $4$ and $2$.