'Poisson model Jacobian of the log-likelihood for each observation
Parameters
params : array-like
The parameters of the model
Returns
score : array-like
The score vector (nobs, k_vars) of the model evaluated at `params`
Notes
.. math:: \frac{\partial\ln L_{i}}{\partial\beta}=\left(y_{i}-\lambda_{i}\right)x_{i}
for observations :math:`i=1,...,n`
where the loglinear model is assumed
.. math:: \ln\lambda_{i}=x_{i}\beta'

def score_obs(self, params):

    offset = getattr(self, 'offset', 0)
    exposure = getattr(self, 'exposure', 0)
    X = self.exog
    L = np.exp(((np.dot(X, params) + offset) + exposure))
    return ((self.endog - L)[:, None] * X)
