'Hessian of NB1 model.'

def _hessian_nb1(self, params):

    if self._transparams:
        alpha = np.exp(params[(-1)])
    else:
        alpha = params[(-1)]
    params = params[:(-1)]
    exog = self.exog
    y = self.endog[:, None]
    mu = self.predict(params)[:, None]
    a1 = (mu / alpha)
    dim = exog.shape[1]
    hess_arr = np.empty(((dim + 1), (dim + 1)))
    dparams = ((exog / alpha) * ((np.log((1 / (alpha + 1))) + special.digamma((y + (mu / alpha)))) - special.digamma((mu / alpha))))
    dmudb = (exog * mu)
    xmu_alpha = ((exog * mu) / alpha)
    trigamma = (special.polygamma(1, ((mu / alpha) + y)) - special.polygamma(1, (mu / alpha)))
    for i in range(dim):
        for j in range(dim):
            if (j > i):
                continue
            hess_arr[(i, j)] = np.sum(((dparams[:, i, None] * dmudb[:, j, None]) + ((xmu_alpha[:, i, None] * xmu_alpha[:, j, None]) * trigamma)), 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) / alpha) * dparams) + (((exog * mu) / alpha) * ((((- trigamma) * mu) / (alpha ** 2)) - (1 / (alpha + 1))))), axis=0)
    hess_arr[(-1), :(-1)] = dldpda
    hess_arr[:(-1), (-1)] = dldpda
    digamma_part = (special.digamma((y + (mu / alpha))) - special.digamma((mu / alpha)))
    log_alpha = np.log((1 / (alpha + 1)))
    alpha3 = (alpha ** 3)
    alpha2 = (alpha ** 2)
    mu2 = (mu ** 2)
    dada = ((((((((((alpha3 * mu) * (((2 * log_alpha) + (2 * digamma_part)) + 3)) - ((2 * alpha3) * y)) + ((alpha2 * mu2) * trigamma)) + (((4 * alpha2) * mu) * (log_alpha + digamma_part))) + (alpha2 * ((2 * mu) - y))) + (((2 * alpha) * mu2) * trigamma)) + (((2 * alpha) * mu) * (log_alpha + digamma_part))) + (mu2 * trigamma)) / ((alpha ** 4) * ((alpha2 + (2 * alpha)) + 1)))
    hess_arr[((-1), (-1))] = dada.sum()
    return hess_arr
