'Returns the score residuals calculated at a given vector of
parameters.
Parameters
params : ndarray
The parameter vector at which the score residuals are
calculated.
Returns
The score residuals, returned as a ndarray having the same
shape as `exog`.
Notes
Observations in a stratum with no observed events have undefined
score residuals, and contain NaN in the returned matrix.'

def score_residuals(self, params):

    surv = self.surv
    score_resid = np.zeros(self.exog.shape, dtype=np.float64)
    mask = np.zeros(self.exog.shape[0], dtype=np.int32)
    w_avg = self.weighted_covariate_averages(params)
    for stx in range(surv.nstrat):
        uft_ix = surv.ufailt_ix[stx]
        exog_s = surv.exog_s[stx]
        nuft = len(uft_ix)
        strat_ix = surv.stratum_rows[stx]
        xp0 = 0.0
        linpred = np.dot(exog_s, params)
        if (surv.offset_s is not None):
            linpred += surv.offset_s[stx]
        linpred -= linpred.max()
        e_linpred = np.exp(linpred)
        at_risk_ix = set([])
        for i in range(nuft)[::(-1)]:
            ix = surv.risk_enter[stx][i]
            at_risk_ix |= set(ix)
            xp0 += e_linpred[ix].sum()
            atr_ix = list(at_risk_ix)
            leverage = (exog_s[atr_ix, :] - w_avg[stx][i, :])
            d = np.zeros(exog_s.shape[0])
            d[uft_ix[i]] = 1
            dchaz = (len(uft_ix[i]) / xp0)
            mrp = (d[atr_ix] - (e_linpred[atr_ix] * dchaz))
            ii = strat_ix[atr_ix]
            score_resid[ii, :] += (leverage * mrp[:, None])
            mask[ii] = 1
            ix = surv.risk_exit[stx][i]
            at_risk_ix -= set(ix)
            xp0 -= e_linpred[ix].sum()
    jj = np.flatnonzero((mask == 0))
    if (len(jj) > 0):
        score_resid[jj, :] = np.nan
    return score_resid
