docs: 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)