MLE for gamma distribution using Bisection method

[GaryBAYLOR]

MLE for gamma distribution using Bisection method

This is a function that I wrote to test bisection method for getting MLE for gamma distribution. After running the function code, we can try an example in R.

> set.seed(101)
> data <- rgamma(500, shape = 4, rate = 2)
> mle.gamma(data)
$hat_alpha
[1] 3.800928

$hat_lambda
[1] 1.88312

We can see that the estimation is reasonable. It can be used for any other estimation when there is no closed form solution. The following is the code for the function mle.gamma().

mle.gamma <- function(X, tol = 1e-6) {
    # X should be a random sample from a gamma distribution
    n <- length(X)
    # define the function f(x) for which we want to find the solution of f(x) = 0
    f <- function(x) n * log(x / mean(X)) + sum(log(X)) - n * digamma(x)

    sign <- function(x) {
        if (x > 0) sign <- 1
        if (x < 0) sign <- -1
        else sign <- 0
        return(sign)
    }

    x.left <- 1
    x.right <- 10
    repeat {
        temp <- (x.left + x.right)/2
        if(abs(f(temp)) < tol) {
            x.final <- temp
            break
        }
        else {
            if(sign(f(temp)) == sign(f(x.right))) x.right <- temp
            else x.left <- temp
        }
    }
    hat_alpha <- x.final
    hat_lambda <- hat_alpha / mean(X)
    lst <- list(hat_alpha = hat_alpha, hat_lambda = hat_lambda)
    return(lst)
}

[GaryBAYLOR]

The MLE by Newton's method

##################################
##  The MLE of gamma distribution
##################################

mle_gamma <- function(x, tol = 1e-4, max_iter = 500) {
    n <- length(x)
    mx <- mean(x)
    g <- function(alpha) {
        log(alpha) - log(mx) + mean(log(x)) - digamma(alpha)
    }
    g_prime <- function(alpha) {
        1 / alpha - trigamma(alpha)
    }

    alpha <- 1
    iter <- 1

    while (iter < max_iter) {
        alpha_new <- alpha - g(alpha) / g_prime(alpha)

        diff <- abs(alpha_new - alpha)
        if(diff < tol | iter > max_iter) break
        alpha <- alpha_new
        iter <- iter + 1
    }

    lambda <- alpha_new / mean(x)   
    res <- c(alpha_new, lambda)
    names(res) <- c("alpha", "lambda")
    print(paste("Iterations: ", iter))
    res
}

> set.seed(101)
> data <- rgamma(500, shape = 4, rate = 2)
> mle.gamma(data)
$hat_alpha
[1] 3.80092760920525

$hat_lambda
[1] 1.88312007132253

> mle_gamma(data)
[1] "Iterations:  7"
           alpha           lambda 
3.80092764816972 1.88312009062697