code: def _score_nbp(y, X, beta, thet, Q):
      r"""
      Negative Binomial Score -- type P likelihood from Greene (2007)
      .. math::
  
          \lambda_i = exp(X\beta)\\
          g_i = \theta \lambda_i^Q \\
          w_i = g_i/(g_i + \lambda_i) \\
          r_i = \theta / (\theta+\lambda_i) \\
          A_i = \left [ \Psi(y_i+g_i) - \Psi(g_i) + ln w_i \right ] \\
          B_i = \left [ g_i (1-w_i) - y_iw_i \right ] \\
          \partial ln \mathcal{L}_i / \partial
              \begin{pmatrix} \lambda_i \\ \theta \\ Q \end{pmatrix}=
              [A_i+B_i]
              \begin{pmatrix} Q/\lambda_i \\ 1/\theta \\ ln(\lambda_i) \end{pmatrix}
              -B_i
              \begin{pmatrix} 1/\lambda_i\\ 0 \\ 0 \end{pmatrix} \\
          \frac{\partial \lambda}{\partial \beta} = \lambda_i \mathbf{x}_i \\
          \frac{\partial \mathcal{L}_i}{\partial \beta} =
              \left (\frac{\partial\mathcal{L}_i}{\partial \lambda_i} \right )
              \frac{\partial \lambda_i}{\partial \beta}
      """
      lamb = np.exp(np.dot(X, beta))
      g = thet * lamb**Q
      w = g / (g + lamb)
      r = thet / (thet + lamb)
      A = digamma(y + g) - digamma(g) + np.log(w)
      B = g * (1 - w) - y * w
      dl = (A + B) * Q / lamb - B * 1 / lamb
      dt = (A + B) * 1 / thet
      dq = (A + B) * np.log(lamb)
      db = X * (dl * lamb)[:, np.newaxis]
      sc = np.array([dt.sum(), dq.sum()])
      sc = np.concatenate([db.sum(axis=0), sc])
      return sc
  
imports: from urllib.request import urlopen
  import numpy
  from numpy.testing import assert_almost_equal
  from scipy.special import digamma
  from scipy.stats import nbinom
  import pandas
  import patsy
  from statsmodels.base.model import GenericLikelihoodModel
  from statsmodels.base.model import GenericLikelihoodModelResults
label_desc: Use image node_id 6 for calling a global function with example usage: _score_nbp(y, X, beta, thet, Q) and returns: sc