code: def approximate_density(
          dist,
          xloc,
          parameters=None,
          cache=None,
          eps=1.e-7
  ):
      
      if parameters is None:
          parameters = dist.prm.copy()
      if cache is None:
          cache = {}
  
      xloc = numpy.asfarray(xloc)
      lo, up = numpy.min(xloc), numpy.max(xloc)
      mu = .5*(lo+up)
      eps = numpy.where(xloc < mu, eps, -eps)*xloc
  
      floc = evaluation.evaluate_forward(
          dist, xloc, parameters=parameters.copy(), cache=cache.copy())
      for d in range(len(dist)):
          xloc[d] += eps[d]
          tmp = evaluation.evaluate_forward(
              dist, xloc, parameters=parameters.copy(), cache=cache.copy())
          floc[d] -= tmp[d]
          xloc[d] -= eps[d]
  
      floc = numpy.abs(floc / eps)
      return floc
docs: Approximate the probability density function.
  
      Args:
          dist : Dist
              Distribution in question. May not be an advanced variable.
          xloc : numpy.ndarray
              Location coordinates. Requires that xloc.shape=(len(dist), K).
          eps : float
              Acceptable error level for the approximations
          retall : bool
              If True return Graph with the next calculation state with the
              approximation.
  
      Returns:
          numpy.ndarray: Local probability density function with
              ``out.shape == xloc.shape``. To calculate actual density function,
              evaluate ``numpy.prod(out, 0)``.
  
      Example:
          >>> distribution = chaospy.Normal(1000, 10)
          >>> xloc = numpy.array([[990, 1000, 1010]])
          >>> print(numpy.around(approximate_density(distribution, xloc), 4))
          [[0.0242 0.0399 0.0242]]
          >>> print(numpy.around(distribution.pdf(xloc), 4))
          [[0.0242 0.0399 0.0242]]