'Jacobian matrix for multinomial logit model log-likelihood
Parameters
params : array
The parameters of the multinomial logit model.
Returns
jac : array-like
The derivative of the loglikelihood for each observation evaluated
at `params` .
Notes
.. math:: \frac{\partial\ln L_{i}}{\partial\beta_{j}}=\left(d_{ij}-\frac{\exp\left(\beta_{j}^{\prime}x_{i}\right)}{\sum_{k=0}^{J}\exp\left(\beta_{k}^{\prime}x_{i}\right)}\right)x_{i}
for :math:`j=1,...,J`, for observations :math:`i=1,...,n`
In the multinomial model the score vector is K x (J-1) but is returned
as a flattened array. The Jacobian has the observations in rows and
the flatteded array of derivatives in columns.'

def score_obs(self, params):

    params = params.reshape(self.K, (-1), order='F')
    firstterm = (self.wendog[:, 1:] - self.cdf(np.dot(self.exog, params))[:, 1:])
    return (firstterm[:, :, None] * self.exog[:, None, :]).reshape(self.exog.shape[0], (-1))
