def quantile(
        self,
        q: Union[float, jnp.ndarray, np.ndarray],
        params: Optional[Params] = None,
        inputs_loader: Optional[InputsLoader] = None,
        mutable: Optional[Mutable] = None,
        calib_params: Optional[CalibParams] = None,
        calib_mutable: Optional[CalibMutable] = None,
        n_target_samples: Optional[int] = 30,
        target_samples: Optional[jnp.ndarray] = None,
        rng: Optional[jax.Array] = None,
        distribute: bool = True,
        **kwargs,
    ) -> Union[float, jnp.ndarray]:
        """
        Estimate the `q`-th quantiles of the likelihood function.

        Parameters
        ----------
        q: Union[float, jnp.ndarray, np.ndarray]
            Quantile or sequence of quantiles to compute, which must be between 0 and 1 inclusive.
        params : Params
            The random parameters of the probabilistic model.
        inputs_loader : InputsLoader
            A loader of input data points.
        mutable : Optional[Mutable]
            The mutable objects used to evaluate the models.
        calib_params : Optional[CalibParams]
            The calibration parameters of the probabilistic model.
        calib_mutable : Optional[CalibMutable]
            The calibration mutable objects used to evaluate the calibrators.
        n_target_samples : int
            Number of target samples to sample for each input data point.
        target_samples: Optional[jnp.ndarray] = None
            Samples of the target variable for each input, used to estimate the quantiles.
        rng: Optional[jax.Array]
            A random number generator. If not passed, this will be taken from the attributes of this class.
        distribute: bool
            Whether to distribute computation over multiple devices, if available.

        Returns
        -------
        jnp.ndarray
            Quantile estimate for each quantile and each input. If multiple quantiles `q` are given, the result's
            first axis is over different quantiles.
        """
        if target_samples is None:
            if params is None or inputs_loader is None:
                raise ValueError(
                    "if `samples` is not passed, then `params` and `inputs_loader` must be passed."
                )
            target_samples = self.sample(
                n_target_samples,
                params,
                inputs_loader,
                mutable,
                calib_params=calib_params,
                calib_mutable=calib_mutable,
                rng=rng,
                distribute=distribute,
                **kwargs,
            )
        return jnp.quantile(target_samples, q, axis=0)

        Estimate the `q`-th quantiles of the likelihood function.

        Parameters
        ----------
        q: Union[float, jnp.ndarray, np.ndarray]
            Quantile or sequence of quantiles to compute, which must be between 0 and 1 inclusive.
        params : Params
            The random parameters of the probabilistic model.
        inputs_loader : InputsLoader
            A loader of input data points.
        mutable : Optional[Mutable]
            The mutable objects used to evaluate the models.
        calib_params : Optional[CalibParams]
            The calibration parameters of the probabilistic model.
        calib_mutable : Optional[CalibMutable]
            The calibration mutable objects used to evaluate the calibrators.
        n_target_samples : int
            Number of target samples to sample for each input data point.
        target_samples: Optional[jnp.ndarray] = None
            Samples of the target variable for each input, used to estimate the quantiles.
        rng: Optional[jax.Array]
            A random number generator. If not passed, this will be taken from the attributes of this class.
        distribute: bool
            Whether to distribute computation over multiple devices, if available.

        Returns
        -------
        jnp.ndarray
            Quantile estimate for each quantile and each input. If multiple quantiles `q` are given, the result's
            first axis is over different quantiles.
        
fortuna/likelihood/regression.py