code: def harrelldavis(x, quantile=0.5, axis=-1):
      """Harrell-Davis robust estimate of the :math:`q^{th}` quantile(s) of the
      data.
  
      .. versionadded:: 0.2.9
  
      Parameters
      ----------
      x : array_like
          Data, must be a one or two-dimensional vector.
      quantile : float or array_like
          Quantile or sequence of quantiles to compute, must be between 0 and 1.
          Default is ``0.5``.
      axis : int
          Axis along which the MAD is computed. Default is the last axis (-1).
          Can be either 0, 1 or -1.
  
      Returns
      -------
      y : float or array_like
          The estimated quantile(s). If ``quantile`` is a single quantile, will
          return a float, otherwise will compute each quantile separately and
          returns an array of floats.
  
      Notes
      -----
      The Harrell-Davis method [1]_ estimates the :math:`q^{th}` quantile by a
      linear combination of  the  order statistics. Results have been tested
      against a Matlab implementation [2]_. Note that this method is also
      used to measure the confidence intervals of the difference between
      quantiles of two groups, as implemented in the shift function [3]_.
  
      See Also
      --------
      plot_shift
  
      References
      ----------
      .. [1] Frank E. Harrell, C. E. Davis, A new distribution-free quantile
         estimator, Biometrika, Volume 69, Issue 3, December 1982, Pages
         635–640, https://doi.org/10.1093/biomet/69.3.635
  
      .. [2] https://github.com/GRousselet/matlab_stats/blob/master/hd.m
  
      .. [3] Rousselet, G. A., Pernet, C. R. and Wilcox, R. R. (2017). Beyond
         differences in means: robust graphical methods to compare two groups
         in neuroscience. Eur J Neurosci, 46: 1738-1748.
         https://doi.org/doi:10.1111/ejn.13610
  
      Examples
      --------
      Estimate the 0.5 quantile (i.e median) of 100 observation picked from a
      normal distribution with zero mean and unit variance.
  
      >>> import numpy as np
      >>> import pingouin as pg
      >>> np.random.seed(123)
      >>> x = np.random.normal(0, 1, 100)
      >>> round(pg.harrelldavis(x, quantile=0.5), 4)
      -0.0499
  
      Several quantiles at once
  
      >>> pg.harrelldavis(x, quantile=[0.25, 0.5, 0.75])
      array([-0.84133224, -0.04991657,  0.95897233])
  
      On the last axis of a 2D vector (default)
  
      >>> np.random.seed(123)
      >>> x = np.random.normal(0, 1, (10, 100))
      >>> pg.harrelldavis(x, quantile=[0.25, 0.5, 0.75])
      array([[-0.84133224, -0.52346777, -0.81801193, -0.74611216, -0.64928321,
              -0.48565262, -0.64332799, -0.8178394 , -0.70058282, -0.73088088],
             [-0.04991657,  0.02932655, -0.08905073, -0.1860034 ,  0.06970415,
               0.15129817,  0.00430958, -0.13784786, -0.08648077, -0.14407123],
             [ 0.95897233,  0.49543002,  0.57712236,  0.48620599,  0.85899005,
               0.7903462 ,  0.76558585,  0.62528436,  0.60421847,  0.52620286]])
  
      On the first axis
  
      >>> pg.harrelldavis(x, quantile=[0.5], axis=0).shape
      (100,)
      """
      x = np.asarray(x)
      assert (
          x.ndim <= 2
      ), "Only 1D or 2D array are supported for this function."
      assert axis in [0, 1, -1], "Axis must be 0, 1 or -1."
  
      # Sort the input array
      x = np.sort(x, axis=axis)
  
      n = x.shape[axis]
      vec = np.arange(n)
      if isinstance(quantile, float):
          quantile = [quantile]
  
      y = []
      for q in quantile:
          # Harrell-Davis estimate of the qth quantile
          m1 = (n + 1) * q
          m2 = (n + 1) * (1 - q)
          w = scipy.stats.beta.cdf(
              (vec + 1) / n, m1, m2
          ) - scipy.stats.beta.cdf((vec) / n, m1, m2)
          if axis != 0:
              y.append((w * x).sum(axis))
          else:
              y.append((w[..., None] * x).sum(axis))  # Store results
  
      if len(y) == 1:
          y = y[
              0
          ]  # Return a float instead of a list if n quantile is 1
      else:
          y = np.array(y)
      return y
  
file_path: pingouin/pingouin/nonparametric.py
imports: import scipy
  import numpy
  import pandas
  from pingouin import remove_na, _check_dataframe, _postprocess_dataframe
label_desc: Use image node_id 8 for calling a global function with example usage: harrelldavis(x, quantile, axis) and returns: y