code: def _stieltjes_approx(dist, order, accuracy, normed, **kws):
      
      kws["rule"] = kws.get("rule", "C")
      assert kws["rule"].upper() != "G"
  
      absisas, weights = chaospy.quad.generate_quadrature(
          accuracy, dist.range(), **kws)
      weights = weights*dist.pdf(absisas)
  
      poly = chaospy.poly.variable(len(dist))
      orth = [poly*0, poly**0]
  
      inner = numpy.sum(absisas*weights, -1)
      norms = [numpy.ones(len(dist)), numpy.ones(len(dist))]
      coeff1 = []
      coeff2 = []
  
      for _ in range(order):
  
          coeff1.append(inner/norms[-1])
          coeff2.append(norms[-1]/norms[-2])
          orth.append((poly-coeff1[-1])*orth[-1] - orth[-2]*coeff2[-1])
  
          raw_nodes = orth[-1](*absisas)**2*weights
          inner = numpy.sum(absisas*raw_nodes, -1)
          norms.append(numpy.sum(raw_nodes, -1))
  
          if normed:
              orth[-1] = orth[-1]/numpy.sqrt(norms[-1])
  
      coeff1.append(inner/norms[-1])
      coeff2.append(norms[-1]/norms[-2])
      coeff1 = numpy.transpose(coeff1)
      coeff2 = numpy.transpose(coeff2)
      norms = numpy.array(norms[1:]).T
      orth = orth[1:]
      return orth, norms, coeff1, coeff2
docs: Stieltjes' method with approximative recurrence coefficients.