<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 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|- - - - - - - - - - .+- - - - . - - - - - - - - - - - | .#| #. | .##| ##. | .##| ##. | .###| ###. M(p) | .###| ###. is the | .###| ###. shaded area | .####| ####. | .####| ####. | ..#####| #####.. | .###########| ###########. t ^ M(p) ^ M(y) | | 1| +++ 1| + | + | + + | ++++++++ | + + | ++ | ++ ++ | ++ | ++ ++ |+++ |+++ +++ 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)