docs: def _estimate_friedrich_coefficients(x, m, r):
      """
      Coefficients of polynomial :math:`h(x)`, which has been fitted to
      the deterministic dynamics of Langevin model
      .. math::
          \dot{x}(t) = h(x(t)) + \mathcal{N}(0,R)
  
      As described by
  
          Friedrich et al. (2000): Physics Letters A 271, p. 217-222
          *Extracting model equations from experimental data*
  
      For short time-series this method is highly dependent on the parameters.
  
      :param x: the time series to calculate the feature of
      :type x: numpy.ndarray
      :param m: order of polynom to fit for estimating fixed points of dynamics
      :type m: int
      :param r: number of quantils to use for averaging
      :type r: float
  
      :return: coefficients of polynomial of deterministic dynamics
      :return type: ndarray
      """
      assert m > 0, "Order of polynomial need to be positive integer, found {}".format(m)
  
      df = pd.DataFrame({'signal': x[:-1], 'delta': np.diff(x)})
      try:
          df['quantiles'] = pd.qcut(df.signal, r)
      except ValueError:
          return [np.NaN] * (m + 1)
  
      quantiles = df.groupby('quantiles')
  
      result = pd.DataFrame({'x_mean': quantiles.signal.mean(), 'y_mean': quantiles.delta.mean()})
      result.dropna(inplace=True)
  
      try:
          return np.polyfit(result.x_mean, result.y_mean, deg=m)
      except (np.linalg.LinAlgError, ValueError):
          return [np.NaN] * (m + 1)