code: def func_quantiles(node, qlist=(.025, .25, .5, .75, .975)):
      
  
      # For very large objects, this will be rather long.
      # Too get the length of the table, use obj.trace.length()
  
      if isinstance(node, Variable):
          func_stacks = node.trace()
      else:
          func_stacks = node
  
      if any(qlist < 0.) or any(qlist > 1.):
          raise TypeError('The elements of qlist must be between 0 and 1')
  
      func_stacks = func_stacks.copy()
  
      N_samp = shape(func_stacks)[0]
      func_len = tuple(shape(func_stacks)[1:])
  
      func_stacks.sort(axis=0)
  
      quants = zeros((len(qlist), func_len), dtype=float)
      alphas = 1. - abs(array(qlist) - .5) / .5
  
      for i in range(len(qlist)):
          quants[i, ] = func_stacks[int(qlist[i] * N_samp), ]
  
      return quants, alphas
docs: Returns an array whose ith row is the q[i]th quantile of the
      function.
  
      :Arguments:
          func_stacks: The samples of the function. func_stacks[i,:]
              gives sample i.
          qlist: A list or array of the quantiles you would like.
  
      :SeeAlso: func_envelopes, func_hist, weightplot