'Probit model score (gradient) vector
Parameters
params : array-like
The parameters of the model
Returns
score : ndarray, 1-D
The score vector of the model, i.e. the first derivative of the
loglikelihood function, evaluated at `params`
Notes
.. math:: \frac{\partial\ln L}{\partial\beta}=\sum_{i=1}^{n}\left[\frac{q_{i}\phi\left(q_{i}x_{i}^{\prime}\beta\right)}{\Phi\left(q_{i}x_{i}^{\prime}\beta\right)}\right]x_{i}
Where :math:`q=2y-1`. This simplification comes from the fact that the
normal distribution is symmetric.'

def score(self, params):

    y = self.endog
    X = self.exog
    XB = np.dot(X, params)
    q = ((2 * y) - 1)
    L = ((q * self.pdf((q * XB))) / np.clip(self.cdf((q * XB)), FLOAT_EPS, (1 - FLOAT_EPS)))
    return np.dot(L, X)
