def _pdf(self, xloc, left, right, cache):
        """
        Probability density function.

        Example:
            >>> print(chaospy.Uniform().pdf([-0.5, 0.5, 1.5, 2.5]))
            [0. 1. 0. 0.]
            >>> print(chaospy.Pow(chaospy.Uniform(), 2).pdf([-0.5, 0.5, 1.5, 2.5]))
            [0.         0.70710678 0.         0.        ]
            >>> print(chaospy.Pow(chaospy.Uniform(1, 2), -1).pdf([0.4, 0.6, 0.8, 1.2]))
            [0.         2.77777778 1.5625     0.        ]
            >>> print(chaospy.Pow(2, chaospy.Uniform()).pdf([-0.5, 0.5, 1.5, 2.5]))
            [0.         0.         0.96179669 0.        ]
            >>> print(chaospy.Pow(2, chaospy.Uniform(-1, 0)).pdf([0.4, 0.6, 0.8, 1.2]))
            [0.         2.40449173 1.8033688  0.        ]
            >>> print(chaospy.Pow(2, 3).pdf([7, 8, 9]))
            [ 0. inf  0.]
        """
        left = evaluation.get_forward_cache(left, cache)
        right = evaluation.get_forward_cache(right, cache)

        if isinstance(left, Dist):
            if isinstance(right, Dist):
                raise StochasticallyDependentError(
                    "under-defined distribution {} or {}".format(left, right))

        elif not isinstance(right, Dist):
            return numpy.inf

        else:

            assert numpy.all(left > 0), "imaginary result"
            x_ = numpy.where(xloc <= 0, -numpy.inf,
                    numpy.log(xloc + 1.*(xloc<=0))/numpy.log(left+1.*(left == 1)))
            num_ = numpy.log(left+1.*(left == 1))*xloc
            num_ = num_ + 1.*(num_==0)

            out = evaluation.evaluate_density(right, x_, cache=cache)/num_
            return out

        x_ = numpy.sign(xloc)*numpy.abs(xloc)**(1./right -1)
        xloc = numpy.sign(xloc)*numpy.abs(xloc)**(1./right)
        pairs = numpy.sign(xloc**right) == 1

        out = evaluation.evaluate_density(left, xloc, cache=cache)
        if numpy.any(pairs):
            out = out + pairs*evaluation.evaluate_density(left, -xloc, cache=cache)

        out = numpy.sign(right)*out * x_ / right
        out[numpy.isnan(out)] = numpy.inf

        return out