code: def quantiles(
      a,
      prob=list([0.25, 0.5, 0.75]),
      alphap=0.4,
      betap=0.4,
      axis=None,
      limit=(),
      masknan=False,
  ):
      """
      Computes empirical quantiles for a data array.
  
      Samples quantile are defined by :math:`Q(p) = (1-g).x[i] +g.x[i+1]`,
      where :math:`x[j]` is the *j*th order statistic, and
      `i = (floor(n*p+m))`, `m=alpha+p*(1-alpha-beta)` and `g = n*p + m - i`.
  
      Typical values of (alpha,beta) are:
          - (0,1)    : *p(k) = k/n* : linear interpolation of cdf (R, type 4)
          - (.5,.5)  : *p(k) = (k+1/2.)/n* : piecewise linear
            function (R, type 5)
          - (0,0)    : *p(k) = k/(n+1)* : (R type 6)
          - (1,1)    : *p(k) = (k-1)/(n-1)*. In this case, p(k) = mode[F(x[k])].
            That's R default (R type 7)
          - (1/3,1/3): *p(k) = (k-1/3)/(n+1/3)*. Then p(k) ~ median[F(x[k])].
            The resulting quantile estimates are approximately median-unbiased
            regardless of the distribution of x. (R type 8)
          - (3/8,3/8): *p(k) = (k-3/8)/(n+1/4)*. Blom.
            The resulting quantile estimates are approximately unbiased
            if x is normally distributed (R type 9)
          - (.4,.4)  : approximately quantile unbiased (Cunnane)
          - (.35,.35): APL, used with PWM ?? JP
          - (0.35, 0.65): PWM   ?? JP  p(k) = (k-0.35)/n
  
      Parameters
      ----------
      a : array_like
          Input data, as a sequence or array of dimension at most 2.
      prob : array_like, optional
          List of quantiles to compute.
      alpha : float, optional
          Plotting positions parameter, default is 0.4.
      beta : float, optional
          Plotting positions parameter, default is 0.4.
      axis : int, optional
          Axis along which to perform the trimming.
          If None (default), the input array is first flattened.
      limit : tuple
          Tuple of (lower, upper) values.
          Values of `a` outside this closed interval are ignored.
  
      Returns
      -------
      quants : MaskedArray
          An array containing the calculated quantiles.
  
      Examples
      --------
      >>> from scipy.stats.mstats import mquantiles
      >>> a = np.array([6., 47., 49., 15., 42., 41., 7., 39., 43., 40., 36.])
      >>> mquantiles(a)
      array([ 19.2,  40. ,  42.8])
  
      Using a 2D array, specifying axis and limit.
  
      >>> data = np.array([[   6.,    7.,    1.],
                           [  47.,   15.,    2.],
                           [  49.,   36.,    3.],
                           [  15.,   39.,    4.],
                           [  42.,   40., -999.],
                           [  41.,   41., -999.],
                           [   7., -999., -999.],
                           [  39., -999., -999.],
                           [  43., -999., -999.],
                           [  40., -999., -999.],
                           [  36., -999., -999.]])
      >>> mquantiles(data, axis=0, limit=(0, 50))
      array([[ 19.2 ,  14.6 ,   1.45],
             [ 40.  ,  37.5 ,   2.5 ],
             [ 42.8 ,  40.05,   3.55]])
  
      >>> data[:, 2] = -999.
      >>> mquantiles(data, axis=0, limit=(0, 50))
      masked_array(data =
       [[19.2 14.6 --]
       [40.0 37.5 --]
       [42.8 40.05 --]],
                   mask =
       [[False False  True]
        [False False  True]
        [False False  True]],
             fill_value = 1e+20)
      """
  
      if isinstance(a, np.ma.MaskedArray):
          return stats.mstats.mquantiles(
              a,
              prob=prob,
              alphap=alphap,
              betap=alphap,
              axis=axis,
              limit=limit,
          )
      if limit:
          marr = stats.mstats.mquantiles(
              a,
              prob=prob,
              alphap=alphap,
              betap=alphap,
              axis=axis,
              limit=limit,
          )
          return ma.filled(marr, fill_value=np.nan)
      if masknan:
          nanmask = np.isnan(a)
          if nanmask.any():
              marr = ma.array(a, mask=nanmask)
              marr = stats.mstats.mquantiles(
                  marr,
                  prob=prob,
                  alphap=alphap,
                  betap=alphap,
                  axis=axis,
                  limit=limit,
              )
              return ma.filled(marr, fill_value=np.nan)
  
      # Initialization & checks ---------
      data = np.asarray(a)
  
      p = np.array(prob, copy=False, ndmin=1)
      m = alphap + p * (1.0 - alphap - betap)
  
      isrolled = False
      # from _quantiles1d
      if axis is None:
          data = data.ravel()  # reshape(-1,1)
          axis = 0
      else:
          axis = np.arange(data.ndim)[axis]
          data = np.rollaxis(data, axis)
          isrolled = True  # keep track, maybe can be removed
  
      x = np.sort(data, axis=0)
      n = x.shape[0]
      returnshape = list(data.shape)
      returnshape[axis] = p
  
      # TODO: check these
      if n == 0:
          return np.empty(len(p), dtype=float)
      elif n == 1:
          return np.resize(x, p.shape)
      aleph = n * p + m
      k = np.floor(aleph.clip(1, n - 1)).astype(int)
      ind = [None] * x.ndim
      ind[0] = slice(None)
      gamma = (aleph - k).clip(0, 1)[ind]
      q = (1.0 - gamma) * x[k - 1] + gamma * x[k]
      if isrolled:
          return np.rollaxis(q, 0, axis + 1)
      else:
          return q
  
file_path: statsmodels/statsmodels/sandbox/stats/stats_mstats_short.py
imports: import numpy
  from numpy import ma
  from scipy import stats
label_desc: Use image node_id 1 for calling a global function with example usage: quantiles(a, prob, alphap, betap, axis, limit, masknan) and returns: stats, ma, np, np, q, ma, np