def _hessian_nb1(self, params):
        if self._transparams: # lnalpha came in during fit
            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

        # for dl/dparams dparams
        dim = exog.shape[1]
        hess_arr = np.empty((dim+1,dim+1))
        #const_arr = a1*mu*(a1+y)/(mu+a1)**2
        # not all of dparams
        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]

        # for dl/dparams dalpha
        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

        # for dl/dalpha dalpha
        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
Determine whether a line is a legal end of RFC 2822 headers.
dataset/ETHPy150Open statsmodels/statsmodels/statsmodels/discrete/discrete_model.py/NegativeBinomial._hessian_nb1--crossover with--dataset/ETHPy150Open ofermend/medicare-demo/socialite/jython/Lib/rfc822.py/Message.islast