'Returns the Hessian of the log partial likelihood evaluated at
`params`, using the Breslow method to handle tied times.'

def breslow_hessian(self, params):

    surv = self.surv
    hess = 0.0
    for stx in range(surv.nstrat):
        uft_ix = surv.ufailt_ix[stx]
        nuft = len(uft_ix)
        exog_s = surv.exog_s[stx]
        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)
        (xp0, xp1, xp2) = (0.0, 0.0, 0.0)
        for i in range(nuft)[::(-1)]:
            ix = surv.risk_enter[stx][i]
            if (len(ix) > 0):
                xp0 += e_linpred[ix].sum()
                v = exog_s[ix, :]
                xp1 += (e_linpred[ix][:, None] * v).sum(0)
                elx = e_linpred[ix]
                xp2 += np.einsum('ij,ik,i->jk', v, v, elx)
            m = len(uft_ix[i])
            hess += (m * ((xp2 / xp0) - (np.outer(xp1, xp1) / (xp0 ** 2))))
            ix = surv.risk_exit[stx][i]
            if (len(ix) > 0):
                xp0 -= e_linpred[ix].sum()
                v = exog_s[ix, :]
                xp1 -= (e_linpred[ix][:, None] * v).sum(0)
                elx = e_linpred[ix]
                xp2 -= np.einsum('ij,ik,i->jk', v, v, elx)
    return (- hess)
