'Hessian of NB2 model.'

def _hessian_nb2(self, params):

    if self._transparams:
        alpha = np.exp(params[(-1)])
    else:
        alpha = params[(-1)]
    a1 = (1 / alpha)
    params = params[:(-1)]
    exog = self.exog
    y = self.endog[:, None]
    mu = self.predict(params)[:, None]
    dim = exog.shape[1]
    hess_arr = np.empty(((dim + 1), (dim + 1)))
    const_arr = (((a1 * mu) * (a1 + y)) / ((mu + a1) ** 2))
    for i in range(dim):
        for j in range(dim):
            if (j > i):
                continue
            hess_arr[(i, j)] = np.sum((((- exog[:, i, None]) * exog[:, j, None]) * const_arr), axis=0)
    tri_idx = np.triu_indices(dim, k=1)
    hess_arr[tri_idx] = hess_arr.T[tri_idx]
    da1 = (- (alpha ** (-2)))
    dldpda = np.sum(((((mu * exog) * (y - mu)) * da1) / ((mu + a1) ** 2)), axis=0)
    hess_arr[(-1), :(-1)] = dldpda
    hess_arr[:(-1), (-1)] = dldpda
    da2 = (2 * (alpha ** (-3)))
    dalpha = (da1 * (((((special.digamma((a1 + y)) - special.digamma(a1)) + np.log(a1)) - np.log((a1 + mu))) - ((a1 + y) / (a1 + mu))) + 1))
    dada = (((da2 * dalpha) / da1) + ((da1 ** 2) * ((((special.polygamma(1, (a1 + y)) - special.polygamma(1, a1)) + (1 / a1)) - (1 / (a1 + mu))) + ((y - mu) / ((mu + a1) ** 2))))).sum()
    hess_arr[((-1), (-1))] = dada
    return hess_arr
