code: def predict_cdf(self, X, quantile, nsamples=200, likelihood_args=()):
          r
          # Get latent function samples
          N = X.shape[0]
          ps = np.empty((N, nsamples))
          fsamples = self._sample_func(X, nsamples)
  
          # Push samples though likelihood cdf
          cdfarg = tuple(chain(atleast_list(self.like_hypers_), likelihood_args))
          for i, f in enumerate(fsamples):
              ps[:, i] = self.likelihood.cdf(quantile, f, *cdfarg)
  
          # Average transformed samples (MC integration)
          p = ps.mean(axis=1)
          p_min = ps.min(axis=1)
          p_max = ps.max(axis=1)
  
          return p, p_min, p_max
docs: r"""
          Predictive cumulative density function of a Bayesian GLM.
  
          Parameters
          ----------
          X : ndarray
              (N*,d) array query input dataset (N* samples, D dimensions).
          quantile : float
              The predictive probability, :math:`p(y^* \leq \text{quantile} |
              \mathbf{x}^*, \mathbf{X}, y)`.
          nsamples : int, optional
              Number of samples for sampling the predictive CDF.
          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*.
          nsamples : int, optional
              The number of samples to draw from the posterior in order to
              approximate the predictive mean and variance.
  
          Returns
          -------
          p : ndarray
             The probability of y* <= quantile for the query inputs, X* of shape
             (N*,).
          p_min : ndarray
              The minimum sampled values of the predicted probability (same shape
              as p)
          p_max : ndarray
              The maximum sampled values of the predicted probability (same shape
              as p)