docs: def calc_finite_diff_terms_for_abc(model_obj,
                                     mle_params,
                                     init_vals,
                                     epsilon,
                                     **fit_kwargs):
      """
      Calculates the terms needed for the finite difference approximations of
      the empirical influence and second order empirical influence functions.
  
      Parameters
      ----------
      model_obj : an instance or sublcass of the MNDC class.
          Should be the model object that corresponds to the model we are
          constructing the bootstrap confidence intervals for.
      mle_params : 1D ndarray.
          Should contain the desired model's maximum likelihood point estimate.
      init_vals : 1D ndarray.
          The initial values used to estimate the desired choice model.
      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.
      fit_kwargs : additional keyword arguments, optional.
          Should contain any additional kwargs used to alter the default behavior
          of `model_obj.fit_mle` and thereby enforce conformity with how the MLE
          was obtained. Will be passed directly to `model_obj.fit_mle`.
  
      Returns
      -------
      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.
  
      References
      ----------
      Efron, Bradley, and Robert J. Tibshirani. An Introduction to the Bootstrap.
          CRC press, 1994. Section 22.6, Equations 22.32 and 22.36.
  
      Notes
      -----
      The returned, symbolic value for `term_minus` does not explicitly appear in
      Equations 22.32 or 22.36. However, it is used to compute a midpoint / slope
      approximation to the finite difference derivative used to define the
      empirical influence function.
      """
      # Determine the number of observations in this dataset.
      num_obs = model_obj.data[model_obj.obs_id_col].unique().size
      # Determine the initial weights per observation.
      init_weights_wide = np.ones(num_obs, dtype=float) / num_obs
      # Initialize wide weights for elements of the second order influence array.
      init_wide_weights_plus = (1 - epsilon) * init_weights_wide
      init_wide_weights_minus = (1 + epsilon) * init_weights_wide
      # Initialize the second order influence array
      term_plus = np.empty((num_obs, init_vals.shape[0]), dtype=float)
      term_minus = np.empty((num_obs, init_vals.shape[0]), dtype=float)
      # Get the rows_to_obs mapping matrix for this model.
      rows_to_obs = model_obj.get_mappings_for_fit()['rows_to_obs']
      # Extract the initial weights from the fit kwargs
      new_fit_kwargs = deepcopy(fit_kwargs)
      if fit_kwargs is not None and 'weights' in fit_kwargs:
          orig_weights = fit_kwargs['weights']
          del new_fit_kwargs['weights']
      else:
          orig_weights = 1
      # Make sure we're just getting the point estimate
      new_fit_kwargs['just_point'] = True
  
      # Populate the second order influence array
      for obs in xrange(num_obs):
          # Note we create the long weights in a for-loop to avoid creating a
          # num_obs by num_obs matrix, which may be a problem for large datasets
          # Get the wide format weights for this observation
          current_wide_weights_plus = init_wide_weights_plus.copy()
          current_wide_weights_plus[obs] += epsilon
  
          current_wide_weights_minus = init_wide_weights_minus.copy()
          current_wide_weights_minus[obs] -= epsilon
          # Get the long format weights for this observation
          long_weights_plus =\
              (create_long_form_weights(model_obj, current_wide_weights_plus,
                                        rows_to_obs=rows_to_obs) * orig_weights)
          long_weights_minus =\
              (create_long_form_weights(model_obj,
                                        current_wide_weights_minus,
                                        rows_to_obs=rows_to_obs) * orig_weights)
          # Get the needed influence estimates.
          term_plus[obs] = model_obj.fit_mle(init_vals,
                                             weights=long_weights_plus,
                                             **new_fit_kwargs)['x']
          term_minus[obs] = model_obj.fit_mle(init_vals,
                                              weights=long_weights_minus,
                                              **new_fit_kwargs)['x']
      return term_plus, term_minus