def compute_pmf(f, x, samples, **kwargs):
    r"""
    Compute the probability mass function given `x` at a range of `x` values
    for :math:`y = f(x|\theta)`

    :math:`P(y|x) = \int P(y=f(x;\theta)|x,\theta) P(\theta) d\theta`

    :math:`\mathrm{pmf}(y|x) = \int_{P(y'|x) < P(y|x)} P(y'|x) dy'`

    Additionally, if a list of log-evidences are passed, along with list of
    functions, samples and optional weights it marginalises over the models
    according to the evidences.

    Parameters
    ----------
    f: function
        function :math:`f(x;\theta)` (or list of functions for each model) with
        dependent variable :math:`x`, parameterised by :math:`\theta`.

    x: 1D array-like
        `x` values to evaluate :math:`f(x;\theta)` at.

    samples: 2D array-like
        :math:`\theta` samples (or list of :math:`\theta` samples) to evaluate
        :math:`f(x;\theta)` at.
        `shape = (nsamples, npars)`

    logZ: 1D array-like, optional
        log-evidences of each model if multiple models are passed.
        Should be same length as the list `f`, and need not be normalised.
        Default: `numpy.ones_like(f)`

    weights: 1D array-like, optional
        sample weights (or list of weights), if desired. Should have length
        same as `samples.shape[0]`.
        Default: `numpy.ones_like(samples)`

    ny: int, optional
        Resolution of y axis.
        Default: `100`

    y: array-like, optional
        Explicit descriptor of `y` values to evaluate.
        Default: `numpy.linspace(min(f), max(f), ny)`

    ntrim: int, optional
        Approximate number of samples to trim down to, if desired. Useful if
        the posterior is dramatically oversampled.
        Default: None

    cache: str, optional
        File root for saving previous calculations for re-use

    parallel, tqdm_args:
        see docstring for :func:`fgivenx.parallel.parallel_apply`


    Returns
    -------
    1D numpy.array:
        `y` values pmf is computed at `shape=(len(y))` or `ny`
    2D numpy.array:
        pmf values at each `x` and `y`  `shape=(len(x),len(y))`
    """

    logZ = kwargs.pop('logZ', None)
    weights = kwargs.pop('weights', None)
    ny = kwargs.pop('ny', 100)
    y = kwargs.pop('y', None)
    ntrim = kwargs.pop('ntrim', 100000)
    parallel = kwargs.pop('parallel', False)
    cache = kwargs.pop('cache', '')
    tqdm_kwargs = kwargs.pop('tqdm_kwargs', {})
    if kwargs:
        raise TypeError('Unexpected **kwargs: %r' % kwargs)

    # y
    if y is not None:
        y = numpy.array(y, dtype='double')
        if len(y.shape) is not 1:
            raise ValueError("y should be a 1D array")

    fsamps = compute_samples(f, x, samples, logZ=logZ,
                             weights=weights, ntrim=ntrim,
                             parallel=parallel, cache=cache,
                             tqdm_kwargs=tqdm_kwargs)

    if y is None:
        ymin = fsamps[~numpy.isnan(fsamps)].min(axis=None)
        ymax = fsamps[~numpy.isnan(fsamps)].max(axis=None)
        y = numpy.linspace(ymin, ymax, ny)

    return y, fgivenx.mass.compute_pmf(fsamps, y, parallel=parallel,
                                       cache=cache, tqdm_kwargs=tqdm_kwargs)
r"""
    Compute the probability mass function given `x` at a range of `x` values
    for :math:`y = f(x|\theta)`

    :math:`P(y|x) = \int P(y=f(x;\theta)|x,\theta) P(\theta) d\theta`

    :math:`\mathrm{pmf}(y|x) = \int_{P(y'|x) < P(y|x)} P(y'|x) dy'`

    Additionally, if a list of log-evidences are passed, along with list of
    functions, samples and optional weights it marginalises over the models
    according to the evidences.

    Parameters
    ----------
    f: function
        function :math:`f(x;\theta)` (or list of functions for each model) with
        dependent variable :math:`x`, parameterised by :math:`\theta`.

    x: 1D array-like
        `x` values to evaluate :math:`f(x;\theta)` at.

    samples: 2D array-like
        :math:`\theta` samples (or list of :math:`\theta` samples) to evaluate
        :math:`f(x;\theta)` at.
        `shape = (nsamples, npars)`

    logZ: 1D array-like, optional
        log-evidences of each model if multiple models are passed.
        Should be same length as the list `f`, and need not be normalised.
        Default: `numpy.ones_like(f)`

    weights: 1D array-like, optional
        sample weights (or list of weights), if desired. Should have length
        same as `samples.shape[0]`.
        Default: `numpy.ones_like(samples)`

    ny: int, optional
        Resolution of y axis.
        Default: `100`

    y: array-like, optional
        Explicit descriptor of `y` values to evaluate.
        Default: `numpy.linspace(min(f), max(f), ny)`

    ntrim: int, optional
        Approximate number of samples to trim down to, if desired. Useful if
        the posterior is dramatically oversampled.
        Default: None

    cache: str, optional
        File root for saving previous calculations for re-use

    parallel, tqdm_args:
        see docstring for :func:`fgivenx.parallel.parallel_apply`


    Returns
    -------
    1D numpy.array:
        `y` values pmf is computed at `shape=(len(y))` or `ny`
    2D numpy.array:
        pmf values at each `x` and `y`  `shape=(len(x),len(y))`