r"""
    Compute the Kullback-Leibler divergence at each value of `x` for the prior
    and posterior defined by `prior_samples` and `samples`.

    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, prior_samples: 2D array-like
        :math:`\theta` samples (or list of :math:`\theta` samples) from
        posterior and prior 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, prior_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)`

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

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

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

    kwargs: further keyword arguments
        Any further keyword arguments are plotting keywords that are passed to
        :func:`fgivenx.plot.plot`.

    Returns
    -------
    1D numpy array:
        dkl values at each value of `x`.

def compute_dkl(f, x, samples, prior_samples, **kwargs):
    r"""
    Compute the Kullback-Leibler divergence at each value of `x` for the prior
    and posterior defined by `prior_samples` and `samples`.

    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, prior_samples: 2D array-like
        :math:`\theta` samples (or list of :math:`\theta` samples) from
        posterior and prior 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, prior_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)`

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

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

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

    kwargs: further keyword arguments
        Any further keyword arguments are plotting keywords that are passed to
        :func:`fgivenx.plot.plot`.

    Returns
    -------
    1D numpy array:
        dkl values at each value of `x`.
    """

    logZ = kwargs.pop('logZ', None)
    weights = kwargs.pop('weights', None)
    prior_weights = kwargs.pop('prior_weights', None)
    ntrim = kwargs.pop('ntrim', None)
    cache = kwargs.pop('cache', '')
    prior_cache = kwargs.pop('prior_cache', '')
    parallel = kwargs.pop('parallel', False)
    tqdm_kwargs = kwargs.pop('tqdm_kwargs', {})

    if kwargs:
        raise TypeError('Unexpected **kwargs: %r' % kwargs)

    if logZ is None:
        logZ = [0]
        f = [f]
        samples = [samples]
        prior_samples = [prior_samples]
        weights = [weights]
        prior_weights = [prior_weights]
        cache = [cache]
        prior_cache = [prior_cache]

    DKLs = []

    for fi, c, pc, s, w, ps, pw in zip(f, cache, prior_cache, samples, weights,
                                       prior_samples, prior_weights):

        fsamps = compute_samples(fi, x, s, weights=w, ntrim=ntrim,
                                 parallel=parallel, cache=c,
                                 tqdm_kwargs=tqdm_kwargs)

        fsamps_prior = compute_samples(fi, x, ps, weights=pw, ntrim=ntrim,
                                       parallel=parallel, cache=pc,
                                       tqdm_kwargs=tqdm_kwargs)

        dkls = fgivenx.dkl.compute_dkl(fsamps, fsamps_prior,
                                       parallel=parallel, cache=c,
                                       tqdm_kwargs=tqdm_kwargs)
        DKLs.append(dkls)

    logZ = numpy.array(logZ)
    DKLs = numpy.array(DKLs)

    Zs = numpy.exp(logZ-logZ.max())
    Zs /= Zs.sum()
    return numpy.sum(Zs * DKLs.transpose(), axis=1)