I am trying to fit a beta distribution to my data using MLE. I calculated the log likelihood function for beta distribution. \begin{equation}\label{eq2} LL(\alpha,\beta) = (\alpha-1)n\bar{\ln x_i}+(\beta-1)n\bar{\ln (1-x_i)}+n\ln \Gamma(\alpha+\beta)-n \ln \Gamma(\alpha)-n \ln \Gamma(\beta) \end{equation} We can ignore $n$ as it is just a scaling factor. In my case values for $\bar{\ln x_i}$ and $\bar{\ln (1-x_i)}$ are -2.62 and -0.15 respectively. So when I plot log likelihood function against the parameter space of $\alpha$ and $\beta$, the function looks concave with a peak around 1 for $\alpha$ and around 5 for $\beta$. So based on these observations, I would conclude that my function is concave for $\alpha$ and $\beta$ around these values.


So far so good. These values also match from what I get by finding the maximum of log likelihood function numerically. But when I calculate the Hessian of this function:
$\begin{bmatrix} \psi_1(\alpha+\beta)-\psi_1(\alpha) & \psi_1(\alpha+\beta)\\ \psi_1(\alpha+\beta) & \psi_1(\alpha+\beta)-\psi_1(\beta) \end{bmatrix}$
where $\psi_1$ is the polygamma function of order 1. It doesn't look like this maxtrix is negative semidefinite. Do you think I am missing something?
Thank you