docs: def quadrect(f, n, a, b, kind='lege', *args, **kwargs):
      """
      Integrate the d-dimensional function f on a rectangle with lower and
      upper bound for dimension i defined by a[i] and b[i], respectively;
      using n[i] points.
  
      Parameters
      ----------
      f : function
          The function to integrate over. This should be a function
          that accepts as its first argument a matrix representing points
          along each dimension (each dimension is a column). Other
          arguments that need to be passed to the function are caught by
          `*args` and `**kwargs`
  
      n : int or array_like(float)
          A length-d iterable of the number of nodes in each dimension
  
      a : scalar or array_like(float)
          A length-d iterable of lower endpoints. If a scalar is given,
          that constant is repeated d times, where d is the number of
          dimensions
  
      b : scalar or array_like(float)
          A length-d iterable of upper endpoints. If a scalar is given,
          that constant is repeated d times, where d is the number of
          dimensions
  
      kind : string, optional(default='lege')
          Specifies which type of integration to perform. Valid
          values are:
  
          lege - Gauss-Legendre
          cheb - Gauss-Chebyshev
          trap - trapezoid rule
          simp - Simpson rule
          N    - Neiderreiter equidistributed sequence
          W    - Weyl equidistributed sequence
          H    - Haber  equidistributed sequence
          R    - Monte Carlo
  
      *args, **kwargs :
          Other arguments passed to the function f
  
      Returns
      -------
      out : scalar (float)
          The value of the integral on the region [a, b]
  
      Notes
      -----
      Based of original function ``quadrect`` in CompEcon toolbox by
      Miranda and Fackler
  
      References
      ----------
      Miranda, Mario J, and Paul L Fackler. Applied Computational
      Economics and Finance, MIT Press, 2002.
  
      """
      if kind.lower() == "lege":
          nodes, weights = qnwlege(n, a, b)
      elif kind.lower() == "cheb":
          nodes, weights = qnwcheb(n, a, b)
      elif kind.lower() == "trap":
          nodes, weights = qnwtrap(n, a, b)
      elif kind.lower() == "simp":
          nodes, weights = qnwsimp(n, a, b)
      else:
          nodes, weights = qnwequi(n, a, b, kind)
  
      out = weights.dot(f(nodes, *args, **kwargs))
      return out