API Reference

def bootstrap(
    data: Any,
    statistic_fn: Union[str, Callable],
    n_boot: int = 1000,
    seed: Optional[int] = None,
    blocksize: Optional[int] = None,
    fn_kwargs: Optional[Dict[str, Any]] = None,
) -> BootstrapDistribution:
    """
    Performs Bayesian bootstrap resampling to estimate a statistic.

    This function performs Bayesian bootstrap resampling by generating `n_boot` resamples from
    the provided `data` and applying the specified statistic function (`statistic_fn`).

    Args:
        data (Any): The data to be resampled. It can be a 1D array, a tuple,
            or a list of arrays where each element represents a different group of data to resample.
        statistic_fn (Union[str, StatisticFunction]): The statistic function to be applied on each
            bootstrap resample. It can either be the name of a registered statistic function or the
            function itself.
        n_boot (int, optional): The number of bootstrap resamples to generate. Default is 1000.
        seed (int, optional): A seed for the random number generator to ensure reproducibility.
            Default is `None`, which means no fixed seed.
        blocksize (int, optional): The block size for resampling. If provided, resampling weights
            are generated in blocks of this size. Defaults to `None`, meaning all resampling weights
            are generated at once.
        fn_kwargs (Dict[str, Any], optional): Additional keyword arguments to be passed to
            the `statistic_fn` for each resample. Default is `None`.

    Returns:
        BootstrapDistribution: An object containing the array with the resampled statistics.

    Raises:
        ValueError: If any data array is not 1D or if the dimensions of the input arrays do not match.

    Example:
        ```python
        data = np.random.randn(100)
        statistic_fn = "mean"
        result = bootstrap(data, statistic_fn)
        print(result)
        print(result.summarize())
        ```

    Notes:
        - The `data` argument can be a single 1D array, or a tuple or list of 1D arrays where each array
          represents a feature of the data.
        - The `statistic_fn` can either be the name of a registered function (as a string) or the function
          itself. If a string is provided, it must match the name of a function in the `statistics.registry`.
        - The function uses the `resample` function to generate bootstrap resamples and apply the statistic
          function to each resample.
    """
    if isinstance(data, np.ndarray):
        if data.ndim != 1:
            raise ValueError(f"Invalid parameter {data.ndim=:}: must be 1.")
        n_data: int = len(data)
    elif isinstance(data, (tuple, list)):
        n_data = len(data[0])
        for i, array in enumerate(data):
            if array.ndim != 1:
                raise ValueError(f"Invalid parameter {data[i].ndim=:}: must be 1.")
            if n_data != len(array):
                raise ValueError(
                    f"Invalid parameter {data[i].shape[0]=:}: must be {n_data=:}."
                )

    if isinstance(statistic_fn, str):
        statistic_fn = get_statistic_fn(statistic_fn)
    estimates = np.array(
        [
            statistic_fn(data=data, weights=weights, **(fn_kwargs or {}))
            for weights in resample(
                n_boot=n_boot,
                n_data=n_data,
                seed=seed,
                blocksize=blocksize,
            )
        ]
    )
    return BootstrapDistribution(estimates=estimates)

def plot(
    bootstrap_distribution: BootstrapDistribution,
    level: float,
    *,
    ax: Optional[plt.Axes] = None,
    n_grid: int = 200,
    label: Optional[str] = None,
    precision: Optional[Union[int, Literal["auto"]]] = None,
) -> plt.Axes:
    """
    Plot the kernel density estimate (KDE) of bootstrap estimates with
    credible interval shading and a vertical line at the mean.

    If an axis is provided, the plot is drawn on it; otherwise, a new figure and axis are created.
    Displays a shaded credible interval and labels the plot with a formatted mean
    and credible interval. If no axis is provided, the figure further is annotated with a title and ylabel,
    ylim[0] positioned at zero, the legend is set, and a tight layout applied.

    Args:
        bootstrap_distribution (BootstrapDistribution): The result of a bootstrap resampling procedure.
        level (float): Credible interval level (e.g., 0.95 for 95% CI).
        ax (plt.Axes, optional): Matplotlib axis to draw the plot on. If None, a new axis is created.
        n_grid (int): Number of grid points to use for evaluating the KDE, default is 200.
        label (str, optional): Optional label for the line. If provided, the label is
            extended to include the mean and credible interval.
        precision (int or "auto" or None, optional): Optional precision for rounding the summary
            values (mean and credible interval). If None (default), no rounding is done; if "auto",
            the precision is computed from the width of the credible interval; if integer, we round to
            this many digits.

    Returns:
        plt.Axes: The axis object containing the plot.
    """
    if ax is None:
        fig, ax = plt.subplots(figsize=(8, 4.5))
    else:
        fig = None

    summary = bootstrap_distribution.summarize(level)

    if precision is not None:
        if precision == "auto":
            summary = summary.round()
        else:
            summary = summary.round(precision)

    param_str = f"{summary.mean} ({summary.ci_low}, {summary.ci_high})"

    if label is not None:
        param_str = f"{label}={param_str}"

    p = gaussian_kde(bootstrap_distribution.estimates)

    x_grid = np.linspace(
        bootstrap_distribution.estimates.min(), bootstrap_distribution.estimates.max(), n_grid
    )
    within_ci = np.logical_and(x_grid >= summary.ci_low, x_grid <= summary.ci_high)
    y_grid = p(x_grid)
    y_mean = p([summary.mean]).item()

    (line,) = ax.plot(x_grid, y_grid, label=param_str)
    color = line.get_color()

    ax.fill_between(
        x_grid[within_ci],
        0,
        y_grid[within_ci],
        facecolor=color,
        alpha=0.5,
    )
    ax.plot([summary.mean] * 2, [0, y_mean], "--", color=color)
    ax.plot([summary.mean], [y_mean], "o", color=color)

    if fig is not None:
        ax.set_title(
            f"Bayesian bootstrap  •  {len(bootstrap_distribution)} resamples, {level * 100:.0f}% CI"
        )
        ax.set_ylim(0, ax.get_ylim()[1])
        ax.set_ylabel("Distribution of estimates")
        ax.legend()
        fig.tight_layout()
    return ax

def get_statistic_fn(name: str) -> StatisticFunction:
    """
    Retrieve a registered statistic function by name.

    Parameters:
        name (str): The lowercase name of the statistic function to retrieve.
            Must be one of:
            - "aggregate"
            - "entropy"
            - "eta_square_dependency"
            - "log_odds"
            - "mean"
            - "median"
            - "mutual_information"
            - "pearson_dependency"
            - "percentile"
            - "probability"
            - "quantile"
            - "self_information"
            - "spearman_dependency"
            - "std"
            - "sum"
            - "variance"

    Returns:
        StatisticFunction: The corresponding function implementation.

    Raises:
        ValueError: If the name does not correspond to a registered function.
    """
    try:
        return STATISTIC_FUNCTIONS[name.lower()]
    except KeyError:
        raise ValueError(
            f"Invalid {name=:}: choose from {list(STATISTIC_FUNCTIONS.keys())}"
        )

def get_statistic_fn_names() -> Tuple[str, ...]:
    """
    Retrieve the names of registered statistic functions.

    Returns:
        Tuple[str, ...]: The names of the available statistic functions.
    """
    return tuple(STATISTIC_FUNCTIONS.keys())

def resample(
    n_boot: int,
    n_data: int,
    seed: Optional[int] = None,
    blocksize: Optional[int] = None,
) -> Generator[FArray, None, None]:
    """
    Generates bootstrap resamples with Dirichlet-distributed weights.

    This function performs resampling by generating weights from a Dirichlet distribution.
    The number of resamples is controlled by the `n_boot` argument, while the size of
    each block of resamples can be adjusted using the `blocksize` argument. The `seed`
    argument allows for reproducible results.

    Args:
        n_boot (int): The total number of bootstrap resamples to generate.
        n_data (int): The number of data points to resample (used for the dimension of the
            Dirichlet distribution).
        seed (int, optional): A random seed for reproducibility (default is `None` for
            random seeding).
        blocksize (int, optional): The number of resamples to generate in each block.
            If `None`, the entire number of resamples is generated in one block.
            Defaults to `None`.

    Yields:
        Generator[FArray, None, None]: A generator that yields each resample
            (a 1D array of floats) as it is generated. Each resample contains Dirichlet-
            distributed weights for the given `n_data`.

    Example:
        ```python
        for r in resample(n_boot=10, n_data=5):
            print(r)
        ```

    Notes:
        - If `blocksize` is specified, the resampling will be performed in smaller blocks,
          which can be useful for parallelizing or limiting memory usage.
        - The function uses NumPy's `default_rng` to generate random numbers, which provides
          a more flexible and efficient interface compared to `np.random.seed`.

    Raises:
        ValueError: If `n_boot` is less than 1 or `n_data` is less than 1.
    """
    if n_boot < 1:
        raise ValueError(f"Invalid parameter {n_boot=:}: must be positive.")
    if n_data < 1:
        raise ValueError(f"Invalid parameter {n_data=:}: must be positive.")
    rng = np.random.default_rng(seed)
    alpha = np.ones(n_data)
    if blocksize is None:
        blocksize = n_boot
    else:
        blocksize = min(max(1, blocksize), n_boot)
    remainder = n_boot
    while remainder > 0:
        size = min(blocksize, remainder)
        weights = rng.dirichlet(alpha=alpha, size=size)
        for w in weights:
            yield w
        remainder -= size

def compute_weighted_aggregate(
    data: FArray,
    weights: FArray,
    *,
    factor: Optional[float] = None,
) -> float:
    """
    Computes a weighted aggregate of the input data.

    This function calculates the dot product of the input `data` and `weights`.
    The function assumes that both `data` and `weights` are 1D arrays of the same
    length and that the weights sum to 1. If a `factor` is provided, the dot product
    is multiplied with it.

    Args:
        data (FArray): A 1D array of numeric values representing the
            data to be aggregated.
        weights (FArray): A 1D array of numeric values representing
            the weights for the data.
        factor (float, optional): A scalar factor to multiply with the computed
            aggregate (default is None).

    Returns:
        float: The computed weighted aggregate, potentially scaled by the `factor`.

    Raises:
        ValueError: If `data` or `weights` are not 1D arrays.
        ValueError: If the shapes of `data` and `weights` do not match.

    Example:
        ```python
        data = np.array([1.0, 2.0, 3.0])
        weights = np.array([0.2, 0.5, 0.3])
        print(compute_weighted_aggregate(data, weights))  # => 2.1
        print(compute_weighted_aggregate(data, weights, factor=1.5))  # => 3.15
        ```

    Notes:
        The weighted aggregate is computed using the dot product between `data` and `weights`.
        The optional `factor` scales the result of this dot product. If no factor is given,
        the aggregation computes the weighted arithmetic mean of the data; if instead the factor
        equals the length of the data array, the aggregation computes the weighted sum.

    """
    validate_array(data=data, weights=weights)
    aggregate = np.dot(weights, data)
    if factor is not None:
        aggregate *= factor
    return aggregate.item()

def compute_weighted_entropy(
    data: IArray,
    weights: FArray,
) -> float:
    """
    Computes a weighted entropy of 1D code data.

    This function calculates the weighted entropy by first computing
    the weighted distribution, dropping the zero elements (contribute
    zero to the following sum), and computing the negative dot product
    between the distribution and log-distribution.

    Args:
        data (IArray): A 1D array of numeric values representing
            the sample data in code format.
        weights (FArray): A 1D array of numeric weights corresponding
            to the data.

    Returns:
        float: The weighted entropy value.

    Raises:
        ValueError: If `data` and `weights` have different shapes or are not 1D.

    Example:
        ```python
        data = np.array([1, 0, 0])
        weights = np.array([0.4, 0.2, 0.4])
        print(compute_weighted_entropy(data, weights))  # => 0.673...
        ```
    """
    validate_array(data=data, weights=weights)
    active_set = get_active_set(data=data)
    distribution = weighted_discrete_distribution(data=active_set, weights=weights)
    distribution = distribution[distribution > 0.0]
    return -np.dot(distribution, np.log(distribution)).item()

def compute_weighted_eta_square_dependency(
    data: IFArray,
    weights: FArray,
) -> float:
    """
    Computes the weighted eta-squared (η²) statistic to assess dependency between
    a categorical and a numerical variable.

    Eta-squared measures the proportion of total variance in the numerical variable
    that is explained by the categorical grouping. It is commonly used in ANOVA-like
    analyses and effect size estimation. The value is bounded between 0 and 1.

    Args:
        data (IFArray): A tuple `(data_cat, data_num)` where:
            - `data_cat` is a 1D array of integer-encoded categorical values.
            - `data_num` is a 1D array of corresponding numeric values.
        weights (FArray): A 1D array of non-negative weights,
            same length as `data_cat`.

    Returns:
        float: Weighted eta-squared value in the range [0, 1], where higher values
        indicate stronger association between the categorical and numeric variable.

    Raises:
        ValueError: If input arrays are not 1D or do not have matching shapes.

    Example:
        ```python
        data_cat = np.array([0, 0, 1, 1])
        data_num = np.array([1.0, 2.0, 3.0, 4.0])
        weights = np.array([0.25, 0.25, 0.25, 0.25)
        print(compute_weighted_eta_square_dependency((data_cat, data_num), weights))  # => 0.8
        ```

    Notes:
        - Internally, η² is computed as the ratio of weighted between-group variance
          to the total weighted variance.
        - The statistic is sensitive to group sizes and imbalance in weights.
        - When all group means equal the global mean, η² is 0.
        - When groups are perfectly separated by the numeric variable, η² is 1.
    """
    data_cat, data_num = data
    mean_sample = compute_weighted_mean(data=data_num, weights=weights)
    group_variance_sum = 0.0
    for value in np.unique(data_cat):
        subset = data_cat == value
        group_weight = weights[subset].sum().item()
        mean_group = compute_weighted_aggregate(
            data=data_num[subset],
            weights=weights[subset],
            factor=1.0 / group_weight,
        )
        group_variance_sum += group_weight * (mean_group - mean_sample) ** 2
    return group_variance_sum / compute_weighted_variance(
        data=data_num,
        weights=weights,
        ddof=0,
    )

def compute_weighted_log_odds(
    data: IArray,
    weights: FArray,
    state: int,
) -> float:
    """
    Computes a weighted log-odds of a state within 1D code data.

    This function calculates the weighted probability of a `state`, `p(state)` by summing
    the `weights` which coincide with `data == state`. The log-odds then
    computes as logarithm of the odds `p(state) / (1 - p(state))`.

    Args:
        data (IArray): A 1D array of numeric values representing
            the sample data in code format.
        weights (FArray): A 1D array of numeric weights corresponding
            to the data.
        state (int): The state for which we estimate the log-odds.

    Returns:
        float: The weighted log-odds value.

    Raises:
        ValueError: If `data` and `weights` have different shapes or are not 1D,
            or if `state` is not in `data`.

    Example:
        ```python
        data = np.array([1, 0, 0])
        weights = np.array([0.4, 0.2, 0.4])
        print(compute_weighted_log_odds(data, weights, state=0))  # => 0.405...
        ```
    """
    probability = compute_weighted_probability(data=data, weights=weights, state=state)
    return math.log(probability / (1.0 - probability))

def compute_weighted_mean(
    data: FArray,
    weights: FArray,
) -> float:
    """
    Computes a weighted mean of the input data.

    This function calculates the weighted arithmetic mean of the input `data`
    and `weights` via the `compute_weighted_aggregate` function.

    Args:
        data (FArray): A 1D array of numeric values representing
            the data to be averaged.
        weights (FArray): A 1D array of numeric values representing
            the weights for the data.

    Returns:
        float: The computed weighted mean.

    Raises:
        ValueError: If `data` or `weights` are not 1D arrays.
        ValueError: If the shapes of `data` and `weights` do not match.

    Example:
        ```python
        data = np.array([1.0, 2.0, 3.0])
        weights = np.array([0.2, 0.5, 0.3])
        print(compute_weighted_mean(data, weights))  # => 2.1
        ```
    """
    return compute_weighted_aggregate(data=data, weights=weights, factor=None)

def compute_weighted_median(
    data: FArray,
    weights: FArray,
    *,
    sorter: Optional[NDArray[np.integer]] = None,
) -> float:
    """
    Computes a weighted median of 1D data using linear interpolation.

    This function calculates the weighted meadian of the given `data` array
    based on the provided `weights` via `compute_weighted_quantile` with parameter
    `quantile=0.5`.

    Args:
        data (FArray): A 1D array of numeric values representing
            the sample data.
        weights (FArray): A 1D array of numeric weights corresponding
            to the data.
        sorter (Optional[NDArray[np.integer]]): Optional array of indices that
            sorts `data`.

    Returns:
        float: The interpolated weighted median value.

    Raises:
        ValueError: If `data` and `weights` have different shapes or are not 1D.

    Example:
        ```python
        data = np.array([1.0, 2.0, 3.0])
        weights = np.array([0.4, 0.2, 0.4])
        print(compute_weighted_median(data, weights))  # => 2.25
        ```
    """
    return compute_weighted_quantile(
        data=data,
        weights=weights,
        quantile=0.5,
        sorter=sorter,
    )

def compute_weighted_mutual_information(
    data: IIArray,
    weights: FArray,
    *,
    normalize: bool = True,
) -> float:
    """
    Computes the weighted mutual information (dependency) between two 1D arrays.

    This function calculates the mutual information dependency between two integer variables
    using the direct mutual information formula on weighted joint distribution estimates.
    The inputs `data_1` and `data_2` are expected to be 1D arrays of the same length, provided
    as a tuple `data`. Each data point is assigned a weight from the `weights` array.

    Args:
        data (IIArray): A tuple of two 1D integer arrays `(data_1, data_2)`
            of equal length.
        weights (FArray): A 1D float array of weights, same length as
            each array in `data`.
        normalize (bool, optional): Whether or not to normalize the mutual information
            to range [0, 1]. Default is True

    Returns:
        float: The weighted mutual information.

    Raises:
        ValueError: If the input arrays are not 1D or have mismatched lengths.

    Example:
        ```python
        data_1 = np.array([0, 0, 1])
        data_2 = np.array([0, 1, 1])
        weights = np.array([0.2, 0.5, 0.3])
        print(compute_weighted_mutual_information((data_1, data_2), weights))  # => 0.146...
        ```

    Notes:
        - The normalized result is bounded by [0, 1], where 0 indicates perfect independence and
          1 indicates perfect dependence.
        - The unnormalized result is bounded between 0 (perfect independence) and
          `min(H(data_0), H(data_1))`, where `H(data_0)` and `H(data_1)` are the entropies of
          the two variables. If we reach the upper bound and `H(data_0) == H(data_1)`,
          we have perfect dependence.
    """
    validate_arrays(data=data, weights=weights)
    active_sets = (
        get_active_set(data[0]),
        get_active_set(data[1]),
    )
    joint_distribution = weighted_discrete_joint_distribution(
        data=active_sets,
        weights=weights,
    )
    distribution_0 = np.sum(joint_distribution, axis=1)
    distribution_1 = np.sum(joint_distribution, axis=0)
    divisor = np.outer(distribution_0, distribution_1).reshape(-1)
    joint_distribution = joint_distribution.reshape(-1)
    non_zero = joint_distribution > 0.0
    mutual_information = np.sum(
        joint_distribution[non_zero]
        * np.log(joint_distribution[non_zero] / divisor[non_zero])
    ).item()
    if mutual_information > 0 and normalize:
        mask = distribution_0 > 0.0
        entropy_0 = -np.dot(distribution_0[mask], np.log(distribution_0[mask])).item()
        mask = distribution_1 > 0.0
        entropy_1 = -np.dot(distribution_1[mask], np.log(distribution_1[mask])).item()
        mutual_information = 2.0 * mutual_information / (entropy_0 + entropy_1)
    return mutual_information

def compute_weighted_pearson_dependency(
    data: FFArray,
    weights: FArray,
    *,
    ddof: int = 0,
) -> float:
    """
    Computes the weighted Pearson correlation coefficient (dependency) between two 1D arrays.

    This function calculates the linear dependency between two variables using a weighted
    version of Pearson's correlation coefficient. The inputs `data_1` and `data_2` are
    expected to be 1D arrays of the same length, provided as a tuple `data`. Each data point
    is assigned a weight from the `weights` array.

    The function normalizes both variables by subtracting their weighted means and dividing
    by their weighted standard deviations, then computes the weighted mean of the element-wise
    product of these normalized arrays.

    Args:
        data (FArray): A tuple of two 1D float arrays `(data_1, data_2)`
            of equal length.
        weights (FArray): A 1D float array of weights, same length as
            each array in `data`.
        ddof (int, optional): Delta degrees of freedom for standard deviation.
            Defaults to 0 (population formula). Use 1 for sample-based correction.

    Returns:
        float: The weighted Pearson correlation coefficient in the range [-1, 1].

    Raises:
        ValueError: If the input arrays are not 1D or have mismatched lengths.

    Example:
        ```python
        data_1 = np.array([1.0, 2.0, 3.0])
        data_2 = np.array([1.0, 2.0, 2.9])
        weights = np.array([0.2, 0.5, 0.3])
        print(compute_weighted_pearson_dependency((data_1, data_2), weights))  # => 0.998...
        ```

    Notes:
        - The function relies on `compute_weighted_mean` and `compute_weighted_std`.
        - The correlation is computed using the formula:
            corr = weighted_mean(z1 * z2)
          where z1 and z2 are the standardized variables.
        - The result is bounded between -1 (perfect negative linear relationship)
          and 1 (perfect positive linear relationship), with 0 indicating no linear dependency.
    """
    data_1, data_2 = data
    weighted_mean_1 = compute_weighted_mean(data=data_1, weights=weights)
    weighted_mean_2 = compute_weighted_mean(data=data_2, weights=weights)
    weighted_std_1 = compute_weighted_std(
        data=data_1,
        weights=weights,
        weighted_mean=weighted_mean_1,
        ddof=ddof,
    )
    weighted_std_2 = compute_weighted_std(
        data=data_2,
        weights=weights,
        weighted_mean=weighted_mean_2,
        ddof=ddof,
    )
    array_1 = (data_1 - weighted_mean_1) / weighted_std_1
    array_2 = (data_2 - weighted_mean_2) / weighted_std_2
    return compute_weighted_mean(data=array_1 * array_2, weights=weights)

def compute_weighted_percentile(
    data: FArray,
    weights: FArray,
    *,
    percentile: float,
    sorter: Optional[NDArray[np.integer]] = None,
) -> float:
    """
    Computes a weighted percentile of 1D data using linear interpolation.

    This function calculates the weighted percentile of the given `data` array
    based on the provided `weights` via `compute_weighted_quantile` with parameter
    `quantile=0.01 * percentile`.

    Args:
        data (FArray): A 1D array of numeric values representing
            the sample data.
        weights (FArray): A 1D array of numeric weights corresponding
            to the data.
        percentile (float): The desired percentile in the interval [0, 100].
        sorter (Optional[NDArray[np.integer]]): Optional array of indices that
            sorts `data`.

    Returns:
        float: The interpolated weighted percentile value.

    Raises:
        ValueError: If `data` and `weights` have different shapes or are not 1D.

    Example:
        ```python
        data = np.array([1.0, 2.0, 3.0])
        weights = np.array([0.2, 0.5, 0.3])
        print(compute_weighted_percentile(data, weights, percentile=70))  # => 2.2
        ```
    """
    return compute_weighted_quantile(
        data=data,
        weights=weights,
        quantile=0.01 * percentile,
        sorter=sorter,
    )

def compute_weighted_probability(
    data: IArray,
    weights: FArray,
    state: int,
) -> float:
    """
    Computes a weighted probability of a state within 1D code data.

    This function calculates the weighted probability of a `state` by summing
    the `weights` which coincide with `data == state`.

    Args:
        data (IArray): A 1D array of numeric values representing
            the sample data in code format.
        weights (FArray): A 1D array of numeric weights corresponding
            to the data.
        state (int): The state for which we estimate the probability.

    Returns:
        float: The weighted probability value.

    Raises:
        ValueError: If `data` and `weights` have different shapes or are not 1D,
            or if `state` is not in `data`.

    Example:
        ```python
        data = np.array([1, 0, 0])
        weights = np.array([0.4, 0.2, 0.4])
        print(compute_weighted_probability(data, weights, state=0))  # => 0.6
        ```
    """
    validate_array(data=data, weights=weights)
    if state not in data:
        raise ValueError(f"Incompatible parameter {state=:}: not included in data.")
    return np.sum(weights[data == state]).item()

def compute_weighted_quantile(
    data: FArray,
    weights: FArray,
    *,
    quantile: float,
    sorter: Optional[IArray] = None,
) -> float:
    """
    Computes a weighted quantile of 1D data using linear interpolation.

    This function calculates the weighted quantile of the given `data` array
    based on the provided `weights`. It uses a normalized cumulative weight
    distribution to determine the interpolated quantile value. The computation
    assumes both `data` and `weights` are 1D arrays of equal length.

    A precomputed `sorter` (array of indices that would sort `data`) can be
    optionally provided to avoid recomputing it internally.

    Args:
        data (FArray): A 1D array of numeric values representing
            the sample data.
        weights (FArray): A 1D array of numeric weights corresponding
            to the data.
        quantile (float): The desired quantile in the interval [0, 1].
        sorter (Optional[NDArray[np.integer]]): Optional array of indices that
            sorts `data`.

    Returns:
        float: The interpolated weighted quantile value.

    Raises:
        ValueError: If `data` and `weights` have different shapes or are not 1D.

    Example:
        ```python
        data = np.array([1.0, 2.0, 3.0])
        weights = np.array([0.2, 0.5, 0.3])
        print(compute_weighted_quantile(data, weights, quantile=0.7))  # => 2.2
        ```

    Notes:
        - If `quantile` is less than or equal to the minimum cumulative weight,
          the smallest data point is returned.
        - If `quantile` is greater than or equal to the maximum cumulative weight,
          the largest data point is returned.
        - Linear interpolation is used between the two closest surrounding data points.
        - Providing a precomputed `sorter` can optimize performance in repeated calls.
    """
    validate_array(data=data, weights=weights)
    if sorter is None:
        sorter = np.argsort(data)
    elif sorter.shape != data.shape:
        raise ValueError(f"Invalid parameter {sorter.shape}: must match {data.shape}.")
    cumulative_weights = np.cumsum(weights[sorter])
    cw_lo, cw_hi = cumulative_weights[[0, -1]]
    cumulative_weights -= cw_lo
    cumulative_weights /= cw_hi - cw_lo
    data = data[sorter]
    if quantile <= cumulative_weights[0]:
        return data[0]
    if quantile >= cumulative_weights[-1]:
        return data[-1]
    idx = np.searchsorted(cumulative_weights, quantile)
    w_lo, w_hi = cumulative_weights[[idx - 1, idx]]
    s_lo, s_hi = data[[idx - 1, idx]]
    w = (quantile - w_lo) / (w_hi - w_lo)
    return (s_lo + w * (s_hi - s_lo)).item()