'Returns a covariance matrix for the proportional hazards model
regresion coefficient estimates that is robust to certain
forms of model misspecification.
Parameters
params : ndarray
The parameter vector at which the covariance matrix is
calculated.
Returns
The robust covariance matrix as a square ndarray.
Notes
This function uses the `groups` argument to determine groups
within which observations may be dependent.  The covariance
matrix is calculated using the Huber-White "sandwich" approach.'

def robust_covariance(self, params):

    if (self.groups is None):
        raise ValueError('`groups`    must    be    specified    to    calculate    the    robust    covariance    matrix')
    hess = self.hessian(params)
    score_obs = self.score_residuals(params)
    grads = {}
    for (i, g) in enumerate(self.groups):
        if (g not in grads):
            grads[g] = 0.0
        grads[g] += score_obs[i, :]
    grads = np.asarray(list(grads.values()))
    mat = grads[None, :, :]
    mat = (mat.T * mat)
    mat = mat.sum(1)
    hess_inv = np.linalg.inv(hess)
    cmat = np.dot(hess_inv, np.dot(mat, hess_inv))
    return cmat
