code: def pdf(self, x_data, step=1e-7):
          
          x_data = numpy.asfarray(x_data)
          shape = x_data.shape
          x_data = x_data.reshape(len(self), -1)
  
          lower, upper = evaluation.evaluate_bound(self, x_data)
          f_data = numpy.zeros(x_data.shape)
          indices = (x_data <= upper) & (x_data >= lower)
          f_data[indices] = evaluation.evaluate_density(self, x_data)[indices]
          f_data = f_data.reshape(shape)
          if len(self) > 1:
              f_data = numpy.prod(f_data, 0)
          return f_data
docs: Probability density function.
  
          If possible the density will be calculated analytically. If not
          possible, it will be approximated by approximating the one-dimensional
          derivative of the forward Rosenblatt transformation and multiplying the
          component parts. Note that even if the distribution is multivariate,
          each component of the Rosenblatt is one-dimensional.
  
          Args:
              x_data (numpy.ndarray):
                  Location for the density function. ``x_data.shape`` must be
                  compatible with distribution shape.
              step (float, numpy.ndarray):
                  If approximation is used, the step length given in the
                  approximation of the derivative. If array provided, elements
                  are used along each axis.
  
          Returns:
              (numpy.ndarray):
                  Evaluated density function values. Shapes are related through
                  the identity ``x_data.shape == dist.shape+out.shape``.