code: def _mom(self, keys, left, right, cache):
          
          if evaluation.get_dependencies(left, right):
              raise evaluation.DependencyError(
                  "sum of dependent distributions not feasible: "
                  "{} and {}".format(left, right)
              )
  
          keys_ = numpy.mgrid[tuple(slice(0, key+1, 1) for key in keys)]
          keys_ = keys_.reshape(len(self), -1)
  
          if isinstance(left, Dist):
              left = [
                  evaluation.evaluate_moment(left, key, cache=cache)
                  for key in keys_.T
              ]
          else:
              left = list(reversed(numpy.array(left).T**keys_.T))
          if isinstance(right, Dist):
              right = [
                  evaluation.evaluate_moment(right, key, cache=cache)
                  for key in keys_.T
              ]
          else:
              right = list(reversed(numpy.array(right).T**keys_.T))
  
          out = numpy.zeros(keys.shape)
          for idx in range(keys_.shape[1]):
              key = keys_.T[idx]
              coef = comb(keys.T, key)
              out += coef*left[idx]*right[idx]*(key <= keys.T)
  
          if len(self) > 1:
              out = numpy.prod(out, 1)
          return out
docs: Statistical moments.
  
          Example:
              >>> print(numpy.around(chaospy.Uniform().mom([0, 1, 2, 3]), 4))
              [1.     0.5    0.3333 0.25  ]
              >>> print(numpy.around(chaospy.Add(chaospy.Uniform(), 2).mom([0, 1, 2, 3]), 4))
              [ 1.      2.5     6.3333 16.25  ]
              >>> print(numpy.around(chaospy.Add(2, chaospy.Uniform()).mom([0, 1, 2, 3]), 4))
              [ 1.      2.5     6.3333 16.25  ]
              >>> print(numpy.around(chaospy.Add(1, 1).mom([0, 1, 2, 3]), 4))
              [1. 2. 4. 8.]