Maximum likelihood parameter estimation in the three-parameter gamma distribution with the use of Mathematica

The number of solutions of the system of the log-likelihood equations for the three-parameter case is still an open problem. Several methods have been developed for finding the solutions. In this article we present a program in Mathematica that can find all the solutions of the system of equations. Furthermore, we examine the case where the global maximum appears at the boundary of the domain of the log-likelihood function and we prove that any consistent estimators appear at the interior with probability tending to one.