docs: def quantile(data, num_breaks):
      """
      Calculate quantile breaks.
  
      Arguments:
      data -- Array of values to classify.
      num_breaks -- Number of breaks to perform.
      """
  
      def scipy_mquantiles(a, prob=list([.25,.5,.75]), alphap=.4, betap=.4, axis=None, limit=()):
          """ function copied from scipy 0.13.3::scipy.stats.mstats.mquantiles """
  
          def _quantiles1D(data,m,p):
              x = numpy.sort(data.compressed())
              n = len(x)
              if n == 0:
                  return numpy.ma.array(numpy.empty(len(p), dtype=float), mask=True)
              elif n == 1:
                  return numpy.ma.array(numpy.resize(x, p.shape), mask=numpy.ma.nomask)
              aleph = (n*p + m)
              k = numpy.floor(aleph.clip(1, n-1)).astype(int)
              gamma = (aleph-k).clip(0,1)
              return (1.-gamma)*x[(k-1).tolist()] + gamma*x[k.tolist()]
  
          # Initialization & checks ---------
          data = numpy.ma.array(a, copy=False)
          if data.ndim > 2:
              raise TypeError("Array should be 2D at most !")
          #
          if limit:
              condition = (limit[0] < data) & (data < limit[1])
              data[~condition.filled(True)] = numpy.ma.masked
          #
          p =  numpy.array(prob, copy=False, ndmin=1)
          m = alphap + p*(1.-alphap-betap)
          # Computes quantiles along axis (or globally)
          if (axis is None):
              return _quantiles1D(data, m, p)
          return numpy.ma.apply_along_axis(_quantiles1D, axis, data, m, p)
  
      return scipy_mquantiles(data, numpy.linspace(1.0 / num_breaks, 1, num_breaks))