docs: def PMF(samples, y):
      """ Compute the probability mass function.
  
          The set of samples defines a probability density P(y),
          which is computed using a kernel density estimator.
  
          From :math:`P(y)` we define:
  
          :math:`\mathrm{pmf}(p) = \int_{P(y)<p} P(y) dy`
  
          This is the cumulative distribution function expressed as a
          function of the probability
  
          We aim to compute :math:`M(y)`, which indicates the amount of
          probability contained outside the iso-probability contour
          passing through :math:`y`::
  
  
               ^ P(y)                   ...
               |                     | .   .
               |                     |       .
              p|- - - - - - - - - - .+- - - - . - - - - - - - - - - -
               |                   .#|        #.
               |                  .##|        ##.
               |                  .##|        ##.
               |                 .###|        ###.     M(p)
               |                 .###|        ###.     is the
               |                 .###|        ###.     shaded area
               |                .####|        ####.
               |                .####|        ####.
               |              ..#####|        #####..
               |          ....#######|        #######....
               |         .###########|        ###########.
               +---------------------+-------------------------------> y
                                    t
  
               ^ M(p)                        ^ M(y)
               |                             |
              1|                +++         1|         +
               |               +             |        + +
               |       ++++++++              |       +   +
               |     ++                      |     ++     ++
               |   ++                        |   ++         ++
               |+++                          |+++             +++
               +---------------------> p     +---------------------> y
              0
  
          Parameters
          ----------
          samples: array-like
              Array of samples from a probability density P(y).
  
          y: array-like (optional)
              Array to evaluate the PDF at
  
          Returns
          -------
          1D numpy.array:
              PMF evaluated at each y value
  
      """
      # Remove any nans from the samples
      samples = numpy.array(samples)
      samples = samples[~numpy.isnan(samples)]
      try:
          # Compute the kernel density estimate
          kernel = scipy.stats.gaussian_kde(samples)
  
          # Add two more samples definitely outside the range and sort them
          mn = min(samples) - 10*numpy.sqrt(kernel.covariance[0, 0])
          mx = max(samples) + 10*numpy.sqrt(kernel.covariance[0, 0])
          y_ = numpy.linspace(mn, mx, len(y)*10)
  
          # Compute the probabilities at each of the extended samples
          ps_ = kernel(y_)
  
          # Compute the masses
          ms = []
          for yi in y:
              # compute the probability at this y value
              p = kernel(yi)
              if p <= max(ps_)*1e-5:
                  m = 0.
              else:
                  # Find out which samples have greater probability than P(y)
                  bools = ps_ > p
  
                  # Compute indices where to start and stop the integration
                  stops = numpy.where(numpy.logical_and(~bools[:-1], bools[1:]))
                  starts = numpy.where(numpy.logical_and(bools[:-1], ~bools[1:]))
  
                  # Compute locations
                  starts = [scipy.optimize.brentq(lambda u: kernel(u)-p,
                                                  y_[i], y_[i+1])
                            for i in starts[0]]
                  starts = [-numpy.inf] + starts
                  stops = [scipy.optimize.brentq(lambda u: kernel(u)-p,
                                                 y_[i], y_[i+1])
                           for i in stops[0]]
                  stops = stops + [numpy.inf]
  
                  # Sum up the masses
                  m = sum(kernel.integrate_box_1d(a, b)
                          for a, b in zip(starts, stops))
              ms.append(m)
          return numpy.array(ms)
  
      except numpy.linalg.LinAlgError:
          return numpy.zeros_like(y)