'Score matrix for multinomial logit model log-likelihood
Parameters
params : array
The parameters of the multinomial logit model.
Returns
score : ndarray, (K * (J-1),)
The 2-d score vector, i.e. the first derivative of the
loglikelihood function, of the multinomial logit model evaluated at
`params`.
Notes
.. math:: \frac{\partial\ln L}{\partial\beta_{j}}=\sum_{i}\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`
In the multinomial model the score matrix is K x J-1 but is returned
as a flattened array to work with the solvers.'

def score(self, params):

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