Predictive percentile interval (upper and lower quantiles).

        Parameters
        ----------
        X : ndarray
            (N*,d) array query input dataset (N* samples, D dimensions).
        percentile : float
            The percentile confidence interval (e.g. 95%) to return.
        nsamples : int, optional
            Number of samples for sampling the predictive percentiles.
        likelihood_args : sequence, optional
            sequence of arguments to pass to the likelihood function. These are
            non-learnable parameters. They can be scalars or arrays of length
            N*.
        multiproc : bool, optional
            Use multiprocessing to paralellise this prediction computation.

        Returns
        -------
        ql : ndarray
            The lower end point of the interval with shape (N*,)
        qu : ndarray
            The upper end point of the interval with shape (N*,)

def predict_interval(self, X, percentile, nsamples=200, likelihood_args=(),
                         multiproc=True):
        """
        Predictive percentile interval (upper and lower quantiles).

        Parameters
        ----------
        X : ndarray
            (N*,d) array query input dataset (N* samples, D dimensions).
        percentile : float
            The percentile confidence interval (e.g. 95%) to return.
        nsamples : int, optional
            Number of samples for sampling the predictive percentiles.
        likelihood_args : sequence, optional
            sequence of arguments to pass to the likelihood function. These are
            non-learnable parameters. They can be scalars or arrays of length
            N*.
        multiproc : bool, optional
            Use multiprocessing to paralellise this prediction computation.

        Returns
        -------
        ql : ndarray
            The lower end point of the interval with shape (N*,)
        qu : ndarray
            The upper end point of the interval with shape (N*,)
        """
        N = X.shape[0]

        # Generate latent function samples per observation (n in N)
        fsamples = self._sample_func(X, nsamples, genaxis=0)

        # Make sure likelihood_args is consistent with work
        if len(likelihood_args) > 0:
            likelihood_args = _reshape_likelihood_args(likelihood_args, N)

        # Now create work for distrbuted workers
        like_hypers = atleast_list(self.like_hypers_)
        work = ((f[0], self.likelihood, like_hypers, f[1:], percentile)
                for f in zip(fsamples, *likelihood_args))

        # Distribute sampling and rootfinding
        if multiproc:
            pool = Pool()
            res = pool.map(_star_rootfinding, work)
            pool.close()
            pool.join()
        else:
            res = [_rootfinding(*w) for w in work]

        # Get results of work
        ql, qu = zip(*res)
        return np.array(ql), np.array(qu)