docs: def calc_2nd_order_influence_abc(mle_params,
                                   term_plus,
                                   term_minus,
                                   epsilon):
      """
      Calculates either a 'positive' finite difference approximation or an
      approximation of a 'positive' finite difference approximation to the the
      2nd order empirical influence array needed to compute the approximate
      boostrap confidence (ABC) intervals. See the 'Notes' section for more
      information on the ambiguous function description.
  
      Parameters
      ----------
      mle_params : 1D ndarray.
          Should contain the desired model's maximum likelihood point estimate.
      term_plus : 2D ndarray.
          Should have one row for each observation. Should have one column for
          each parameter in the parameter vector being estimated. Elements should
          denote the finite difference term that comes from adding a small value
          to the observation corresponding to that elements respective row.
      term_minus : 2D ndarray.
          Should have one row for each observation. Should have one column for
          each parameter in the parameter vector being estimated. Elements should
          denote the finite difference term that comes from subtracting a small
          value to the observation corresponding to that elements respective row.
      epsilon : positive float.
          Should denote the 'very small' value being used to calculate the
          desired finite difference approximations to the various influence
          functions. Should be 'close' to zero.
  
      Returns
      -------
      second_order_influence : 2D ndarray.
          Should have one row for each observation. Should have one column for
          each parameter in the parameter vector being estimated. Elements should
          denote the second order empirical influence of the associated
          observation on the associated parameter.
  
      References
      ----------
      Efron, Bradley, and Robert J. Tibshirani. An Introduction to the Bootstrap.
          CRC press, 1994. Section 22.6, Equation 22.36.
  
      Notes
      -----
      This function is based on the code in Efron's original Bootstrap library
      written in S-plus. It is not equivalent to the 'positive' finite difference
      approximation to Equation 22.36, where epsilon is set to a small float. The
      reason for this discrepancy is becaue `term_minus` is not equal to the
      third term in the numerator in Equation 22.36. That term uses
      `(1 - epsilon)P^0` whereas `term_minus` uses `(1 + epsilon)P^0`. At the
      limit, both of these terms would be `P^0` and therefore equal. I think
      Efron's original code was making the assumption that the two terms are
      approximately equal to conserve computational resources. Either that or
      Equation 22.36, as printed, is incorrect because its third term really
      should be `(1 + epsilon)P^0`.
      """
      # Calculate the denominator of each row of the second order influence array
      denominator = epsilon**2
      # Initialize the second term used to calculate each row in the second order
      # influence array
      term_2 = np.broadcast_to(2 * mle_params, term_plus.shape)
      # Calculate the second order empirical influence array.
      second_order_influence = np.zeros(term_plus.shape, dtype=float)
      # Only perform the calculations if the values in the component terms are
      # sufficiently different. This is to prevent floating point errors. atol is
      # simply set to a 'reasonably small' number.
      diff_idx = ~np.isclose(term_plus + term_minus, term_2, atol=1e-12, rtol=0)
      if diff_idx.any():
          second_order_influence[diff_idx] =\
              ((term_plus[diff_idx] - term_2[diff_idx] + term_minus[diff_idx]) /
               denominator)
      return second_order_influence