'For computing marginal effects standard errors.
This is used only in the case of discrete and count regressors to
get the variance-covariance of the marginal effects. It returns
[d F / d params] where F is the predicted probabilities for each
choice. dFdparams is of shape nobs x (J*K) x (J-1)*K.
The zero derivatives for the base category are not included.
Transform can be \'dydx\' or \'eydx\'. Checking is done in margeff
computations for appropriate transform.'

def _derivative_predict(self, params, exog=None, transform='dydx'):

    if (exog is None):
        exog = self.exog
    if (params.ndim == 1):
        params = params.reshape(self.K, (self.J - 1), order='F')
    eXB = np.exp(np.dot(exog, params))
    sum_eXB = (1 + eXB.sum(1))[:, None]
    (J, K) = lmap(int, [self.J, self.K])
    repeat_eXB = np.repeat(eXB, J, axis=1)
    X = np.tile(exog, (J - 1))
    F0 = (((- repeat_eXB) * X) / (sum_eXB ** 2))
    F1 = (((eXB.T[:, :, None] * X) * (sum_eXB - repeat_eXB)) / (sum_eXB ** 2))
    F1 = F1.transpose((1, 0, 2))
    other_idx = (~ np.kron(np.eye((J - 1)), np.ones(K)).astype(bool))
    F1[:, other_idx] = ((((- eXB.T[:, :, None]) * X) * repeat_eXB) / (sum_eXB ** 2)).transpose((1, 0, 2))[:, other_idx]
    dFdX = np.concatenate((F0[:, None, :], F1), axis=1)
    if ('ey' in transform):
        dFdX /= self.predict(params, exog)[:, :, None]
    return dFdX
