docs: def calc_quadratic_coef_abc(model_object,
                              mle_params,
                              init_vals,
                              empirical_influence,
                              std_error,
                              epsilon,
                              **fit_kwargs):
      """
      Calculates the quadratic coefficient needed to compute the approximate
      boostrap confidence (ABC) intervals.
  
      Parameters
      ----------
      model_object : 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.
      empirical_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 empirical influence of the associated observation on the
          associated parameter.
      std_error : 1D ndarray.
          Contains the standard error of the MLE estimates for use in the ABC
          confidence intervals.
      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
      -------
      quadratic_coef : 1D ndarray.
          Contains a measure of nonlinearity of the MLE estimation function as
          one moves in the 'least favorable direction.'
  
      References
      ----------
      Efron, Bradley, and Robert J. Tibshirani. An Introduction to the Bootstrap.
          CRC press, 1994. Section 22.6, Equation 22.37.
      """
      # Determine the number of observations in this dataset.
      num_obs = float(empirical_influence.shape[0])
      # Calculate the standardized version of the empirical influence values
      # assuming mean zero. This is the mean of empirical influence values
      # (sum over n) divided by the standard error of the mean of the empiricial
      # influence values, i.e. the square root of [n^2 * variance] = n * std_err
      # where std_err is the standard error of the sampling distribution of the
      # influence values.
      standardized_influence =\
          empirical_influence / (num_obs**2 * std_error[None, :])
      # Determine the initial weights per observation.
      init_weights_wide = (np.ones(int(num_obs), dtype=float) / num_obs)[:, None]
      # Create the wide weights for the various terms of the quadratic_coef
      term_1_wide_weights =\
          (1 - epsilon) * init_weights_wide + epsilon * standardized_influence
      term_3_wide_weights =\
          (1 - epsilon) * init_weights_wide - epsilon * standardized_influence
      # Get the rows_to_obs mapping matrix for this model.
      rows_to_obs = model_object.get_mappings_for_fit()['rows_to_obs']
      # Initialize the terms for the quadratic coefficients.
      expected_term_shape = (init_vals.shape[0], init_vals.shape[0])
      term_1_array = np.empty(expected_term_shape, dtype=float)
      term_3_array = np.empty(expected_term_shape, dtype=float)
      # 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
      # Calculate the various terms of the quadratic_coef
      for param_id in xrange(expected_term_shape[0]):
          # Get a 'long-format' array of the weights per observation
          term_1_long_weights =\
              (create_long_form_weights(model_object,
                                        term_1_wide_weights[:, param_id],
                                        rows_to_obs=rows_to_obs) * orig_weights)
          term_3_long_weights =\
              (create_long_form_weights(model_object,
                                        term_3_wide_weights[:, param_id],
                                        rows_to_obs=rows_to_obs) * orig_weights)
          # Populate the given row of the term_1 and term_3 arrays.
          term_1_array[param_id] =\
              model_object.fit_mle(init_vals,
                                   weights=term_1_long_weights,
                                   **new_fit_kwargs)['x']
          term_3_array[param_id] =\
              model_object.fit_mle(init_vals,
                                   weights=term_3_long_weights,
                                   **new_fit_kwargs)['x']
      # Extract the desired terms from their 2d arrays
      term_1 = np.diag(term_1_array)
      term_3 = np.diag(term_3_array)
      term_2 = 2 * mle_params
      # Calculate the quadratic coefficient
      quadratic_coef = np.zeros(term_1.shape, dtype=float)
      denominator = epsilon**2
      # Only perform the calculations if the values in the component terms are
      # sufficiently different. This is to prevent floating point errors. Note
      # that atol is simply set to a 'reasonably small' number.
      diff_idx = ~np.isclose(term_1 + term_3, term_2, atol=1e-10, rtol=0)
      if diff_idx.any():
          quadratic_coef[diff_idx] =\
              ((term_1[diff_idx] - term_2[diff_idx] + term_3[diff_idx]) /
               denominator)
      return quadratic_coef