import mpmath


def nll_invhess(x, k, theta):
    """
    Gamma distribution inverse of the hessian of the negative log-likelihood.
    """
    k = mpmath.mpf(k)
    theta = mpmath.mpf(theta)

    N = len(x)
    sumx = mpmath.fsum(x)
    # sumlnx = mpmath.fsum(mpmath.log(t) for t in x)

    dk2 = -N*mpmath.psi(1, k)
    dkdtheta = -N/theta
    dtheta2 = -2*sumx/theta**3 + N*k/theta**2

    det = dk2*dtheta2 - dkdtheta**2

    return mpmath.matrix([[-dtheta2/det, dkdtheta/det],
                          [dkdtheta/det, -dk2/det]])