def prediction_interval_coverage_probability(
    lower_bounds: Array, upper_bounds: Array, targets: Array
) -> Array:
    """
    Compute the prediction interval coverage probability (PICP). This is the fraction of data points for which the
    true targets lie within the estimated interval. This is supported only for scalar target data.

    Parameters
    ----------
    lower_bounds: Array
        Predictive lower bounds. These are the lower bounds of the estimated predictive intervals. This can either
        be a one-dimensional array with entry corresponding to different data points, or a two-dimensional array
        with first axis corresponding to different data points, and second axis with only one dimension.
    upper_bounds: Array
        Predictive upper bounds. These are the upper bounds of the estimated predictive intervals. This can either
        be a one-dimensional array with entry corresponding to different data points, or a two-dimensional array
        with first axis corresponding to different data points, and second axis with only one dimension.
    targets: Array
        A two-dimensional array of target variables, or a one-dimensional array with second dimension of size 1.

    Returns
    -------
    Array
        The computed PICP.
    """
    if targets.ndim == 1:
        targets = targets[:, None]
    elif targets.shape[1] > 1:
        raise ValueError(
            """This metric is supported only for target data such that `target.shape[1] == 1`, but
        `target.shape[1] == {}` was found.""".format(
                targets.shape[1]
            )
        )
    if upper_bounds.ndim == 1:
        upper_bounds = upper_bounds[:, None]
    if upper_bounds.shape[1] != 1:
        raise ValueError(
            """The second axis of `upper_bounds` must contain only one component."""
        )
    if lower_bounds.ndim == 1:
        lower_bounds = lower_bounds[:, None]
    if lower_bounds.shape[1] != 1:
        raise ValueError(
            """The second axis of `lower_bounds` must contain only one component."""
        )

    satisfies_upper = targets <= upper_bounds
    satisfies_lower = targets >= lower_bounds
    return jnp.mean(satisfies_lower * satisfies_upper)

    Compute the prediction interval coverage probability (PICP). This is the fraction of data points for which the
    true targets lie within the estimated interval. This is supported only for scalar target data.

    Parameters
    ----------
    lower_bounds: Array
        Predictive lower bounds. These are the lower bounds of the estimated predictive intervals. This can either
        be a one-dimensional array with entry corresponding to different data points, or a two-dimensional array
        with first axis corresponding to different data points, and second axis with only one dimension.
    upper_bounds: Array
        Predictive upper bounds. These are the upper bounds of the estimated predictive intervals. This can either
        be a one-dimensional array with entry corresponding to different data points, or a two-dimensional array
        with first axis corresponding to different data points, and second axis with only one dimension.
    targets: Array
        A two-dimensional array of target variables, or a one-dimensional array with second dimension of size 1.

    Returns
    -------
    Array
        The computed PICP.
    
fortuna/metric/regression.py