docs: def _find_statistics(X, y, variogram_function,
                       variogram_model_parameters, coordinates_type):
      """Calculates variogram fit statistics.
      Returns the delta, sigma, and epsilon values for the variogram fit.
      These arrays are used for statistics calculations.
  
      Parameters
      ----------
      X: ndarray
          float array [n_samples, n_dim], the input array of coordinates
      y: ndarray
          float array [n_samples], the input array of measurement values
      variogram_function: callable
          function that will be called to evaluate variogram model
      variogram_model_parameters: list
          user-specified parameters for variogram model
      coordinates_type: str
          type of coordinates in X array, can be 'euclidean' for standard
          rectangular coordinates or 'geographic' if the coordinates are lat/lon
  
      Returns
      -------
      delta: ndarray
          residuals between observed values and kriged estimates for those values
      sigma: ndarray
          mean error in kriging estimates
      epsilon: ndarray
          residuals normalized by their mean error
      """
  
      delta = np.zeros(y.shape)
      sigma = np.zeros(y.shape)
  
      for i in range(y.shape[0]):
  
          # skip the first value in the kriging problem
          if i == 0:
              continue
  
          else:
              k, ss = _krige(X[:i, :], y[:i], X[i, :], variogram_function,
                             variogram_model_parameters, coordinates_type)
  
              # if the estimation error is zero, it's probably because
              # the evaluation point X[i, :] is really close to one of the
              # kriging system points in X[:i, :]...
              # in the case of zero estimation error, the results are not stored
              if np.absolute(ss) < eps:
                  continue
  
              delta[i] = y[i] - k
              sigma[i] = np.sqrt(ss)
  
      # only use non-zero entries in these arrays... sigma is used to pull out
      # non-zero entries in both cases because it is guaranteed to be positive,
      # whereas delta can be either positive or negative
      delta = delta[sigma > eps]
      sigma = sigma[sigma > eps]
      epsilon = delta/sigma
  
      return delta, sigma, epsilon