'Score vector for NB2 model'

def _score_nbin(self, params, Q=0):

    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 = ((1 / alpha) * (mu ** Q))
    if Q:
        dparams = (((exog * mu) / alpha) * ((np.log((1 / (alpha + 1))) + special.digamma((y + (mu / alpha)))) - special.digamma((mu / alpha))))
        dalpha = (((alpha * ((y - (mu * np.log((1 / (alpha + 1))))) - (mu * ((special.digamma((y + (mu / alpha))) - special.digamma((mu / alpha))) + 1)))) - (mu * ((np.log((1 / (alpha + 1))) + special.digamma((y + (mu / alpha)))) - special.digamma((mu / alpha))))) / ((alpha ** 2) * (alpha + 1))).sum()
    else:
        dparams = (((exog * a1) * (y - mu)) / (mu + a1))
        da1 = (- (alpha ** (-2)))
        dalpha = ((((((special.digamma((a1 + y)) - special.digamma(a1)) + np.log(a1)) - np.log((a1 + mu))) - ((a1 + y) / (a1 + mu))) + 1).sum() * da1)
    if self._transparams:
        return np.r_[(dparams.sum(0), (dalpha * alpha))]
    else:
        return np.r_[(dparams.sum(0), dalpha)]
