def _get_efron_values_single(self, X, T, E, weights, beta):
        """
        Calculates the first and second order vector differentials, with respect to beta.
        Note that X, T, E are assumed to be sorted on T!

        A good explanation for Efron. Consider three of five subjects who fail at the time.
        As it is not known a priori that who is the first to fail, so one-third of
        (φ1 + φ2 + φ3) is adjusted from sum_j^{5} φj after one fails. Similarly two-third
        of (φ1 + φ2 + φ3) is adjusted after first two individuals fail, etc.

        From https://cran.r-project.org/web/packages/survival/survival.pdf:

        "Setting all weights to 2 for instance will give the same coefficient estimate but halve the variance. When
        the Efron approximation for ties (default) is employed replication of the data will not give exactly the same coefficients as the
        weights option, and in this case the weighted fit is arguably the correct one."

        Parameters
        ----------
        X: array
            (n,d) numpy array of observations.
        T: array
            (n) numpy array representing observed durations.
        E: array
            (n) numpy array representing death events.
        weights: array
            (n) an array representing weights per observation.
        beta: array
            (1, d) numpy array of coefficients.

        Returns
        -------
        hessian:
            (d, d) numpy array,
        gradient:
            (1, d) numpy array
        log_likelihood: float
        """

        n, d = X.shape
        hessian = np.zeros((d, d))
        gradient = np.zeros((d,))
        log_lik = 0

        # Init risk and tie sums to zero
        x_death_sum = np.zeros((d,))
        risk_phi, tie_phi = 0, 0
        risk_phi_x, tie_phi_x = np.zeros((d,)), np.zeros((d,))
        risk_phi_x_x, tie_phi_x_x = np.zeros((d, d)), np.zeros((d, d))

        # Init number of ties and weights
        weight_count = 0.0
        tied_death_counts = 0
        scores = weights * np.exp(np.dot(X, beta))

        # Iterate backwards to utilize recursive relationship
        for i in range(n - 1, -1, -1):
            # Doing it like this to preserve shape
            ti = T[i]
            ei = E[i]
            xi = X[i]
            score = scores[i]
            w = weights[i]

            # Calculate phi values
            phi_i = score
            phi_x_i = phi_i * xi
            phi_x_x_i = np.outer(xi, phi_x_i)

            # Calculate sums of Risk set
            risk_phi = risk_phi + phi_i
            risk_phi_x = risk_phi_x + phi_x_i
            risk_phi_x_x = risk_phi_x_x + phi_x_x_i

            # Calculate sums of Ties, if this is an event
            if ei:
                x_death_sum = x_death_sum + w * xi
                tie_phi = tie_phi + phi_i
                tie_phi_x = tie_phi_x + phi_x_i
                tie_phi_x_x = tie_phi_x_x + phi_x_x_i

                # Keep track of count
                tied_death_counts += 1
                weight_count += w

            if i > 0 and T[i - 1] == ti:
                # There are more ties/members of the risk set
                continue
            elif tied_death_counts == 0:
                # Only censored with current time, move on
                continue

            # There was atleast one event and no more ties remain. Time to sum.
            #
            # This code is near identical to the _batch algorithm below. In fact, see _batch for comments.
            #
            weighted_average = weight_count / tied_death_counts

            if tied_death_counts > 1:
                increasing_proportion = np.arange(tied_death_counts) / tied_death_counts
                denom = 1.0 / (risk_phi - increasing_proportion * tie_phi)
                numer = risk_phi_x - np.outer(increasing_proportion, tie_phi_x)
                a1 = np.einsum("ab,i->ab", risk_phi_x_x, denom) - np.einsum(
                    "ab,i->ab", tie_phi_x_x, increasing_proportion * denom
                )
            else:
                denom = 1.0 / np.array([risk_phi])
                numer = risk_phi_x
                a1 = risk_phi_x_x * denom

            summand = numer * denom[:, None]
            a2 = summand.T.dot(summand)

            gradient = gradient + x_death_sum - weighted_average * summand.sum(0)

            log_lik = log_lik + np.dot(x_death_sum, beta) + weighted_average * np.log(denom).sum()
            hessian = hessian + weighted_average * (a2 - a1)

            # reset tie values
            tied_death_counts = 0
            weight_count = 0.0
            x_death_sum = np.zeros((d,))
            tie_phi = 0
            tie_phi_x = np.zeros((d,))
            tie_phi_x_x = np.zeros((d, d))

        return hessian, gradient, log_lik
Calculates the first and second order vector differentials, with respect to beta.
        Note that X, T, E are assumed to be sorted on T!

        A good explanation for Efron. Consider three of five subjects who fail at the time.
        As it is not known a priori that who is the first to fail, so one-third of
        (φ1 + φ2 + φ3) is adjusted from sum_j^{5} φj after one fails. Similarly two-third
        of (φ1 + φ2 + φ3) is adjusted after first two individuals fail, etc.

        From https://cran.r-project.org/web/packages/survival/survival.pdf:

        "Setting all weights to 2 for instance will give the same coefficient estimate but halve the variance. When
        the Efron approximation for ties (default) is employed replication of the data will not give exactly the same coefficients as the
        weights option, and in this case the weighted fit is arguably the correct one."

        Parameters
        ----------
        X: array
            (n,d) numpy array of observations.
        T: array
            (n) numpy array representing observed durations.
        E: array
            (n) numpy array representing death events.
        weights: array
            (n) an array representing weights per observation.
        beta: array
            (1, d) numpy array of coefficients.

        Returns
        -------
        hessian:
            (d, d) numpy array,
        gradient:
            (1, d) numpy array
        log_likelihood: float